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Introduction 



Quite a few physical problems give rise to quantum lattice models. Among these are 
descriptions of high-Tc superconducting materials, metals, insulators and magnets. 
A quantum lattice model is characterised by quantum objects — spins, fermions or 
bosons — nicely positioned on a regular lattice. We will restrict ourselves to spins. 
These are located at the lattice points. The low temperature behaviour of such a 
model is strongly influenced by quantum zero point fluctuations. These usually sup- 
press the classical order. Sometimes the classical order is even destroyed giving room 
for other, non-classical types of ordering. 

In this thesis insight in the low temperature behaviour will be obtained by a 
numerical study of the ground state at zero temperature. Two schemes are possible 
for this purpose and both will be employed. 

The first scheme makes a direct estimate of the ground state wave function. An 
approximate ground state is built and the properties are analysed afterwards. The 
most successful member of this class is the Density Matrix Renormalisation Group 
(DMRG). It consists of a systematic, iterative procedure. At every iteration step 
the energy is minimised in a given subspace of the configurational Hilbert space. 
This gives a variational approximation to the ground state. Afterwards, a part of the 
subspace is enlarged and another part is truncated. This transformation is tuned to 
keep an optimal fraction of the approximation within the altered subspace. To preview 
one of the comparisons with the second scheme, the DMRG does not suffer from the 
'sign-problem' that hampers many other approaches. Below the sign-problem will be 
explained. 

White introduced the DMRG in 1992 and it has proven to be extremely suc- 
cessful for one-dimensional quantum models [11^, R3, Bl] and two-dimensional classi- 



cal models 0, |T^. Applications to higher dimensional quantum models are relatively 
T3| , ^ ^ 0. In this thesis we investigate what can be achieved in two 



rare 



dimensions and find that the method has substantial more difficulty with the two- 
dimensional geometry. We will try to explain this limitation. 

The second scheme does not attempt to approximate the ground state but aims at 
a direct sampling of the properties instead. From this class we will employ the Green 
Function Monte Carlo simulation (GFMC) [Tl|, Dimensionality does not 



play such an important role here as it does in the DMRG method. The essential 
assumption of the Green Function Monte Carlo simulation, like that of any Monte 
Carlo simulation, is that the properties of a system can be obtained by measuring 
them in many representative configurations. Every measurement Xa is accompanied 
by a weight Mq, to express its importance and the average over these measurements, 
J2a XaMal ^a, will yield the properties. Green Function Monte Carlo simulations 
suffer from two important shortcomings. The most important one is that quantum 
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mechanical models that contain either frustration or fermions require an extension of 
these weights M to negative values and the average Y.a ^aMa/ Y.a becomes prone 
to noise as the individual, positive and negative weights are exponentially larger 
that their sum. This complication is called the sign-problem. The second limitation 
is that a priori a good estimate of the ground state has to exist to help the simulation 
distinguish relevant configurations from less relevant ones. With the introduction 
of the Fixed Node Monte Carlo (FNMC) the sign-problem can also be cured by 
incorporating a very good approximation to the ground state in the simulation. The 
extension of the Fixed Node Monte Carlo to the Green Function Monte Carlo with 
Stochastic Reconfiguration (GFMCSR) |^ makes the end result less dependent on 
the quality of the approximation. Still a good approximation remains vital for the 
method. 

The objective of this thesis is two-fold. First, we study the DMRG method to 
understand its capabilities and limitations in two dimensions. Second, the DMRG 
will be integrated with a Green Function Monte Carlo simulation. DMRG provides a 
systematic approximation to the ground state based on the energy. The correlation 
functions are biased by the implementation of the method. This shortcoming can 
be overcome by a Green Function Monte Carlo simulation using the DMRG wave 
function as a guiding wave function. This thesis introduces this new and promising 
combination of the DMRG and a Green Function Monte Carlo simulation for the first 
time. 

To achieve these objectives, we will first study a well-known model, the Ising 
model in a Transverse Field (ITF) to analyse the DMRG method. Afterwards the 
DMRG is applied to an unsolved problem, the two-dimensional, antiferromagnetic, 
frustrated Heisenberg model. Although this model has been attacked by a variety 
of methods |]14|, |T6|, ^ 0, ^ |5^ no definite results exist so far for the 
quantum phase diagram of this model as Monte Carlo simulations are hampered 
by a fundamental problem called the sign-problem. Strong indications for the exis- 
tence of a phase without classical ordering are found by DMRG calculations and the 
combined effort of the DMRG and the Green Function Monte Carlo with Stochastic 
Reconfiguration. The next paragraphs describe this brief outline in more detail. 

The first two chapters focus on the method itself. Chapter one introduces the 
two-dimensional Ising model in a Transverse Field. It has a direct and clear mapping 
to a highly anisotropic three-dimensional Ising model, making it almost a blue-print 
for a model with a quantum phase transition. It does not suffer from the sign-problem 
and cluster Monte Carlo simulations have yielded high quality numerical results |^ to 
which we can compare our results. DMRG can only handle relatively small systems 
and finite-size scaling techniques are necessary to extend the results to larger system 
sizes. A large fraction of the first chapter is devoted to developing these finite-size 
scaling techniques. Thanks to the power of the DMRG combined with these scaling 
relations we can numerically establish the critical field of the two-dimensional Ising 
model in a Transverse Field upto three significant figures. 

The main subject of chapter two is the DMRG technique. A new variant of 
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the method is introduced and afterwards the general properties are described. It 
seems that DMRG will need amendments or modifications for larger two-dimensional 
systems. In combination with scaling relations we can however establish the two- 
dimensional behaviour. 

In the following two chapters the DMRG technique is apphed to the frustrated 
Heisenberg model. The frustration in that model appears by competition of nearest- 
neighbour and next nearest-neighbour interactions, an the consequences for the phase 
diagram are unclear. The essential question is whether a phase with no classical equiv- 
alent exists for intermediate range of frustration. Chapter three is dedicated to the 
description of the model and the introduction of the spin stiffness ps- The spin stiffness 
is an excellent indicator of long-range order and should reveal whether an intermediate 
phase exists. It is studied using a mean field approximation in the Schwinger-boson 
representation. The reason to resort to a relatively complex mean field approximation 
is that the ground state of the frustrated Heisenberg model is rotationally invariant. 
The mean fields can thus not be simply the spin expectation values ( {S) = 0). The 
Schwinger-Boson Mean Field approximation (SBMF) inserts the mean fields in the 
interactions of neighbouring spins, Si ■ Sj. This not only overcomes the complication 
of the rotational invariance, but it also extends the mean fields to local correlations. 
It yields a rotationally invariant ground state but the spin length ^ cannot be strictly 
conserved. 

The spin stiffness in this approximation does not reveal an intermediate phase, 
but it serves very well for an analysis of the finite-size scaling behaviour. This helps 
us to extend the numerical results of the next chapter to an infinitely large system. 

Chapter four outlines the numerical calculation including the technique to obtain 
the spin stiffness, guided by the results of chapter three. Like many other methods. 
The DMRG do not give a definitive answer on the existence of an intermediate phase, 
but it provides clear information on infinitely long strips of widths upto eight sites. 

The final chapter combines the DMRG and the Green Function Monte Carlo with 
Stochastic Reconfiguration. The frustrated Heisenberg model belongs to the class 
where Monte Carlo simulations suffer from the sign-problem. As mentioned above 
a good guidance is essential for the Green Function Monte Carlo with Stochastic 
Reconfiguration. For a long time finding an proper approximation to the ground 
state has been the bottleneck of all Green Function Monte Carlo variants as a large 
amount of research time had to be spent on designing it. The DMRG can provide 
such an approximate ground state for many different models, including the frustrated 
Heisenberg model. In this chapter it is outlined how to combine both methods and 
the phase diagram of a 10 x 10 system is studied. This combination of the DMRG and 
a Green Function Monte Carlo simulation is new and promising. Further extensions, 
along the lines of forward- walking schemes |Q, may even be able to obtain accurate 
values for the spin-spin correlations. 

This thesis will hopefully provide a good understanding of the intricacies of the 
DMRG method. The last chapter resolves a long standing problem of the Green 
Function Monte Carlo and the last three chapters give indications of the intermedi- 
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ate phase of the frustrated Heisenberg model although no definite statements can be 
made. 
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1 Ising model in a Transverse field 



1.1 Introduction 

Since the beginning of the 1960s, the Ising model in a Transverse field (ITF) has been 
studied. In first instance, this quantum mechanical model was employed to describe 
the order-disorder transition in some double- well ferroelectric systems like KH2P04^ 
crystals. This interest has survived to the present day, but the scope has widened. 

A decade later the renormalisation group and with it the notion of universality 
was introduced. The Ising model often served as a test ground and consequently it 
was scrutinised. The d dimensional ITF can be mapped onto a d + 1 dimensional 
Ising model. This relation makes it an excellent vehicle to introduce the concepts of 
phase transitions in the realm of quantum mechanics. 



Sachdev, Read and others |^2|, |13[ have used the ITF for the same role as the Ising 
model has played in the context of classical critical phenomena: a blue-print of phase 
transitions and universality. Maybe the most important difference lies in the fact that 
it is not the temperature T that induces a phase transition, but a coupling constant 
H that can drastically alter the properties of the system. With the disappearance 
of the temperature, T = 0, it is the ground state that exhibits this quantum phase 
transition. 

On the outset our intention is to investigate the density matrix renormalisation 
group (DMRG). The ITF is chosen as a 'toy-model' both because of its rich behaviour 
and its simple description. 

In this chapter, the model will be explained. The exact results in one dimension 
(a chain) will be reviewed and subsequently a large effort is made to uncover the 
finite-size scaling behaviour for the two-dimensional case. Given the restriction of the 
DMRG, which will be discussed in the next chapter, it is worthwhile to scale the 
length of the system to infinity first, after which the finite width can be scaled away. 
The results, table |1-3|, clearly support such a two-stage process. 



A good review of the ITF has been published by Chakrabarti et al. |T3[, relieving 
us of the duty to go into great detail. The numerical treatment is left to the next 
chapter. 



1.2 The Model 

Consider a two dimensional square lattice with length L and width W. The lattice is 
periodic in both directions and each lattice site contains a spin-^. The Hamiltonian 
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is given by 

L W 

niTF = E E {-^^U^r+i,.. + + '^Hsi^} (1-2.1) 

1=1 w=l 

where the iSf^ are the usual spin-1 matrices satisfying the commutation relations 

["S^n^^ ^i/,w'] = i^i,v^w,w'eap^Sl^ , a, /3, 7 = X, z (1.2.2) 

Clearly the energy scale of this Hamiltonian is set to unity. We have chosen the 
convention of spin-1, in contrast with the Pauli- matrices frequently used in this field. 
The factor of 4 and 2 are inserted in the definition (|1.2.1| ) to make our results directly 
comparable to the main part of the literature. 

Let us summarise the main symmetry properties of this model. The Hamiltonian 
with field H can be transformed into one with —H by a spin rotation round the 
(S'^-axis over 180 degrees. This is a unitary transformation so we have the freedom to 
choose H > 0. The model is translation and reflection symmetric in both directions. 

An important symmetry that will be extensively used, is the spin-reversal operator 
S. It is associated with a rotation over 180 degrees round the axis; 

s^s'^s = -s'' , s^s^s = -sy , s^s's = s\ 

and can be expressed as 

5 = exp ^7r|^E'5z% + ^j • (1.2.3) 

The offset of LW/2 allows us to associate the quantum number S = 1 with the 
ground state \iPq) for different system sizes. ( In short: S\iPq) = iV'o)). One can state 
that the spin reversal operator samples the number of up-spins and returns whether 
it is even, 5 = 1, or odd, S = —1. 

If if —> 0, we end up with a simple, classical, two-dimensional Ising model. The 
ground state is degenerate; all spins point forwards or backwards in the iS^-direction. 
The associated phase is named the classically ordered phase. The two classical so- 
lutions, forward and backward pointing spins, can be superposed in a quantum me- 
chanical sense. In this manner states can be obtained that are either even [S = +1) 
in spin- reversal or odd (5 = —1). The statement made before about the ground state, 
S\ipo) = \ipo) still holds for one of them, but the other ground state {ipi) lies in the 
odd subspace, S\ipi) = — iV'i). 

In the other extreme, H — ^ oo, the degeneracy is lifted. The model essentially 
describes free spins in an external field. The ground state is unique and has all 
spins pointing down in the iS^-direction. This is the reference state for the quantum 
disordered phase and again has the quantum number 5 = +1. The lowest excitation 
differs from the ground state by the reversal of one spin. So it belongs to the class 
iS = — 1. We will extensively study the energy gap A between the lowest excitation 
(in S = —1) and the ground state (in S = +1). A = Ei — Eq. 



1.3 Connection with the 3-dimensional Ising model 
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Figure 1-1. The phase diagram for the ITF and the energy gap A associated with it. At 
the critical field He a phase transition occurs and the gap A opens up. Further explanation 
is given in the text. 



There is a phase transition between the classically ordered and the quantum 
disordered state. A clear signature of this phase transition is the appearance of the 



gap A, which occurs for a critical value H = H^. In figure |1-1| all these properties are 
summarised in a graphical representation. As we will show later the relation between 
the gap A and the field can be made more explicit by 



\H-Hr 



;i.2.4) 



for H > He- Below the critical field, H < the ground state becomes degenerate. 
The gap A should then be defined as the energy difference between the ground state 
and the first excitation within the even subspace for equation ( |1.2.4| ) to hold. In the 
rest of the chapter we will not redefine the gap but only consider H > He. 

It will take another section to prove that the critical exponent u is identical to 
the critical exponent of the three-dimensional classical Ising model. After that we 
focus on extracting both the critical field and the critical exponent u for the 
two-dimensional ITF. 



1.3 Connection with the 3-dimensional Ising model 

The ITF has an intimate connection with a highly anisotropic Ising model in one 
dimension higher. Although we will only show this explicitly for the 2-dimensional 
ITF, the argument holds in general. 
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As starting point we take a highly anisotropic, ferromagnetic Ising model in three 
dimensions. The classical spins = ±1 interact according to the Hamiltonian 

- [3n,las = K^Y. sis] + ^11 E sis], (1.3.1) 

on a cubic lattice at an inverse temperature The first term contains the interac- 
tions between nearest neighbours (denoted by [ij]) in the plane perpendicular to the 
anisotropy axis. The second term contains those along the anisotropy direction. The 
coupling constants and K\\ give rise to this anisotropy. We link them to the field 



in equation (|1.2.1| ) where the energy scale was set to unity, 

K^ = e , e-2^ii =eE , e < 1, (1.3.2) 
is a small, positive parameter and K\\ is a large one. 

The partition function can be calculated by the transfer matrix method. We choose 
the transfer direction along the anisotropy axis. The Hamiltonian l-idas can be split 
in interactions between slices perpendicular to the anisotropy direction. For each pair 
of adjacent slices we then obtain a local, transfer Hamiltonian Titrans for all slices. 
The connections with the ITF will be found in the partition sum. 



N 

e 

n=l 



The connection to the quantum mechanical world is made by introducing summations 
over complete bases |q;„) for the slices in this product; 



N 



^ = E nK|e-^'^"-"">„+i), 

|o„> n=l 

with periodic boundary conditions, |aAr+i) = \ai). The parameters i^j, and K\\ were 
chosen so anisotropic that we may approximate the exponentials; 



n\ exp I K\\ E sl-s] = e^ii(a„|l + 2e '^^^^Y.^ijf^n+i). 

exp I E s!s] \an+i) = («n|l + 4:K^Y.^l,wiSl+l,w + '5z%+i)|«n+l) 
/ l,w 



All of the above is combined in the resulting partition function 

N f T L W 

Ze-^^ii = Tr n 1 + ^ E E (45,%(5f_,i,^ + 5,%^^) + 2HSl^ 

ri=l \ 11=1 ■w=l 

The right side of the equation can easily be mapped onto the definition of the ITF 
in equation ( p..2.1| ) by some simple, unitary transformations: we exchange the S^- 
and (S^-direction and afterwards rotate the spins over vr round the axis (using the 



1.4 Exact solutions for the ITF chain 
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unitary operator exp(27rX]j ^j))- The trace is invariant under these unitary transfor- 
mation, so the outcome is 

Ze-^^w = Tie-^''^'^^. (1.3.3) 

Careful analysis reveals that this relation holds up to order 0{e'^). Along with this 
connection come quite a few others; the correlation function ^, formally defined by 
the covariance of Ising spins located at location (/, w, 1) and (Z, w, n); 

can be related to the gap A in the ITF by inserting complete bases in the same 
fashion as above. The result is 

= In ^ = e{E, - Eo) = eA. (1.3.4) 



As mentioned before in the paragraph after equation (|1.2.4|) this relation only holds 



for H > He- Below the critical field the gap has to be redefined. 

The reduced temperature t = {T — Tc)/Tc in the classical model is replaced by the 
reduced field h = [H — He) / He- The relation for critical exponent u thus also finds 
its equivalent; 

- ~ ^ A ~ h'r 

In this expression we use that for the 3D Ising model z = 1 



1.4 Exact solutions for the ITF chain 

The properties of the two-dimensional system can only be obtained by use of finite- 
size scaling. In doing so, the results for the open ITF chain will provide a crucial 
reference. Here we briefiy review those. The reason not to follow Pfeuty ||3^ is that 



he considers an infinitely long chain, whereas here it is essential that the chain is both 
finite and open. Given the beauty of exact solutions, we can not resist in presenting 
also the results for the finite periodic chain. 

The essential ingredient to the exact solution of the ITF chain is a Jordan- Wigner 
transformation |]5S[ to spin- less fermions Cj, 



5+ = ne^^^^c,. 
In terms of these operators the Hamiltonian reads 

L L-l 

n = -HL + c]c, -JY. (c] - C,) (c]+i + C,+i) + H^er- (1.4. 1) 

i=i j=i 
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The term "Hper governs the interactions between the first and the last site. If the chain 
is open, these do not exist, TCper = 0. For a periodic chain 

Hper = j{ci- Cl) {cl + Cq) S. (1.4.2) 

This equation reintroduces another observable, the spin-reversal operator S, defined 
as 

L 

S = exp in ^ c]cj. 



The current definition is identical to the original one in equation ( |1.2.3| ). It is a 



conserved quantity as the only change in the number of fermions occurs in the second 
term of the Hamiltonian ( |1.4.1|) by pair creation or annihilation. As an empty chain 
corresponds to all spins pointing downwards, the ground state has to be in the even, 
S = 1, subspace. Moreover these pairs of operators in the Hamiltonian Ti make it 
already clear that it can be diagonalised via a Bogoliubov transformation. For the 
periodic chain we can go even further and derive exact expressions for the excitation 
spectrum. The open-chain properties require the numerical diagonalisation of a L x L 
matrix to get all energy levels and eigenstates. 

1.4.1 The periodic chain 

We want to Fourier transform the particle creation and annihilation operators in the 
usual manner 

L 
1=1 

The allowed fc-values depend on the number of fermions present on the chain; for an 
odd number , S = —1, the boundary term TCper, equation (|1.4.2| ), has the same sign 
as all other interactions. The Hamiltonian Ti. becomes translational invariant and we 
can take the regular fc-values; 

A; = ^ , 0<l <L. (1.4.3) 

For even number of fermions, 5=1, the boundary conditions are antiperiodic and 
the set of fc- values has to be adjusted accordingly; 

k = E±Dl,0<l<L. (1.4.4) 

The Hamiltonian Ti. is now converted to the momentum space representation and to 
represent it compactly we define p{k) = —2H + 2Jcos(/c) and q{k) = 2Jsin(/c). It 
reads 



n =-HN - pik){clck + cUc. 



-k 

0<k<Tr 

J J 



0<fc<7r 



1.4 Exact solutions for the ITF chain 
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Length 


S=-l 


S=l 


even 




A; = 0,7r 


odd 


k = n 


k = 



Table 1-1. The k = and k = tt values will only occur for certain length and number of 
fermions. . 



where the allowed fc- values are defined in ( |1.4.3| ) and (|1.4.4| ). The wave vector k = 
and k = n — treated in ?io,7r — are not always allowed. If they are, see table |1-1| , they 
play a special role in that they already appear in a diagonal form; 

^0,. = ^0 allowed P(0)coCo + 5^ allowed P(^)4c^- 

Apart from these, all the terms have to undergo a Bogoliubov transformation. Next we 
provide the main formulas to diagonalise those: define the canonical transformation 

Cfc = u{k)r]k - iv{k)r]lf., 

with 

u{k) = cos 9k , v{k) = sinOk, 
cos 2^, = ^ , A(A;) = ^ (k) + (k) . 

The final expression for the Hamiltonian now is 

0<fc<7r 

with 

k 

It is important to stress that the form of the Hamiltonian Ti. depends on the number 
of fermions present. Creating one excitation by creating a single fermion will change 
the parity of the number of particles; the only allowed excitations within the subspace 
are built from pairs of fermions, i. e. ''7|'7]|V'o)- To obtain an excitation in the other 
subspace, we have to start with the lowest energy state within that subspace and 
create pairs of excitations on that. 



1.4.2 The open chain 

In the case of the periodic chain we are fortunate that the translational invariance 
allows a Fourier transformation. The modes become almost decoupled afterwards and 
an exact expression can be derived by a Bogoliubov transformation. 

The open chain is not translation invariant and no further analytical steps can be 
taken. The numerical procedure remains straightforward. Van Hemmen p4[ describes 
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15 



10 



Open Ising chain 

Closed Ising chain 




Figure 1-2. The behaviour of the energy gap A as function of the combination of x = H—Hc 
and the length L of the system. The only difference between the lines is the boundary 
conditions. 



in detail how to derive both the energy spectrum and all eigenstates of a fermionic 
system for a general pair- wise interaction. Applying the relevant transformations leads 
again to a diagonal Hamiltonian containing excitation occupation numbers offsetted 
by the ground state energy Eq. 

In figure |1-2| the scaled gap LA, on which we will focus from now on, is depicted 
for both the open and the closed chain. They show similar behaviour, although in 
the vicinity of the phase transition their ratio becomes quite large. 



1.5 Finite-size scaling 

We want to establish properties of the 2D ITF through numerical means. In general it 
is not possible to obtain direct quantitative information on an infinitely large system. 
A widely used indirect route is to calculate the desired quantities for a set of finite 
systems and extract their values for the infinite-size limit by investigating their de- 
pendencies on the systems size. This method is called finite-size scaling. Anticipating 
that for our numerical method, the density matrix renormalisation group, the length 
of the system imposes less restrictions than the width, we will first develop the scaling 
relations to scale the length L —>■ oo. Afterwards we derive the relations for W —>■ oo. 
A final subsection is spent on studying the scaling behaviour for fixed aspect ratio 
L/W and L — »• cxo. The properties of the 2D system are independent of the route 
taken to derive them, therefore both approaches should yield the same results in the 
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limit of infinite system size. 

Although the Density Matrix Renormalisation Group (DMRG) will not be intro- 
duced until the next chapter, the results are incorporated in the following scaling 
analysis. DMRG calculations were performed to obtain the lowest energies in the 
even and odd subspace. The difference is the energy gap A. The parameters of the 
calculations were H = 3.00, . . . , 3.10 in steps of 0.01 and L/W = 2, . . . , 5 in steps of 
0.5. For widths W = 7,8 the largest ratio was L/W = 4. 

1.5.1 1-Dimensional, W fixed, L ^ oo 

There are two reasons to scale the length L to infinity first. The first one has to do 
with the open boundary conditions. As we will discuss in the next chapter, closing 
the boundaries in the length direction and making the system translation invari- 
ant will severely hamper the accuracy of the calculation. Open systems give rise to 
complex finite-size effects by their lack of translational invariance. These effects will 
disappear if we scale the length L oo where open and periodic systems become 
indistinguishable. 

The second reason is of a more practical nature. The fact that DMRG functions 
so extremely well for quantum chains |51[], indicates that varying the length L of the 
system will not have a large impact on the accuracy. Unfortunately this does not 
hold for the width W of the system. In the next chapter, section (|2.6.1D , we will show 
that in order to maintain the accuracy, the size of the calculation grows exponentially 
with the width W of the system. With a large range of lengths and only a few widths 
that we can handle, it seems wise to remove all the dependence on the length first. 
We will thereby obtain fairly accurate results for infinitely long strips. 

Let us now outline the procedure: once the length L has become sufficiently large, 
a system of dimensions L x W, with fixed width W, will start to behave as a one- 
dimensional system. Such a system will by arguments of universality resemble an 
open chain, so we expect for fixed width W, that 

LA{x, W, L) = A{W) i^h{B{W)xL) + ^h{B{W)xL) + . . .| . (1.5.1) 

This scaling relation contains quite some new elements that deserve either an intro- 
duction or an explanation: 

• LA. Given the connection to the Ising model, A has to correspond to an inverse 
length. Rescaling the length will thus automatically lead to the appearance of 
the combination LA. 

• /o(. . .) and /i(...)/L. These are respectively the leading order and the first 
correction of the finite-size behaviour of the open ITF chains. We obtained 
these functions numerically. Knowing that the correction arises from a surface 
contribution (the chain is open), the factor 1/L in /i(. . .)/L is obvious. 
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Figure 1-3. The critical field Hc(W) can be obtained by 'smoothing' the curve of LA versus 
xL. Depicted is the case for W = 5 and Hc(W) = 2.867. 



• xL. Universal behaviour concerns the system properties in the vicinity of the 
phase transition. Every width has its own critical field Hc{W), so x = H — 
Hc{W). One dimensional critical behaviour is accompanied by the critical ex- 
ponent u = 1. The combination thus has to be xL'^ = xL. 

• A{W) and B{W). The theory of universality makes predictions on the exponents 
near a phase transition, not on the overall scales. These can vary for the different 
width strips. 

It has to be stressed that only the form of the leading term, fo{B{W)xL), is fixed 
by universality. The form of the correction term, fi{B{W)xL)/L, is an assumption. 
We just take the simplest form that also suits the open chain. 

The scales v4(iy), B{W) and the critical fields Hc(W) can be considered fitting 
parameters to make relation (|1.5.1| ) agree with the data. In the following we will 
describe how we obtained these parameters and what the outcome is. 

The standard fitting procedure is to adjust the parameters such that all data 
points are fitted best according to a least-square functional. The critical field Hc{W) 
can also be derived independently through another, simpler route. In practice this 
approach is taken. If the right Hc{W) is chosen, plotting LA versus xL must give 
a smooth curve for large enough L. From this feature we will extract Hc{W). An 
example of this procedure is given in figure |1-3| . 

With Hc{W) already found, only y4(iy) and B{W) have to be obtained from data 
collapse. We can fit the data x, L, A to the formula ( p..5.1|) , as both /o and /i can be 
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Figure 1-4. The data obtained by DMRG calculations is fitted to scaling relation (1.5.1) 
The figure shows scaling collapse nicely. 



obtained from the open ITF chain. Figure |l-4] reveals that indeed the data shows nice, 
universal behaviour with the calculated /q. The corrections due to /i are removed 



from this figure. They are relatively small. Table L-2 lists the results. 

The general expression for the gap A at infinite system length L = oo can be 
found by examining the behaviour of /o(|/) and fi{y) at large argument y more 
closely: for large argument y (or field) the ITF Hamiltonian describes free spins in an 
external magnetic field with strength y. An excitation can be made by fiipping one 
spin upwards, thus 

lim = 1, 

y^oo y 

with no corrections dependent on size. Therefore the gap A will become 

A{x,W)= lim A{x,W,L) = A{W)B{W)[H - H^{W)]. (1.5.2) 

L— >oo 

This relation for the gap A of a finite-width strip and the numerical value of the 
quantities appearing in it, table are the results of this subsection. In the light of 
the gap expression, the trends in the numerics become clearer: 

The critical field Hc{W) is monotonically increasing and should approach the 
two-dimen-sional value, 

lim H,{W) = H,. 
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w 


H,{W) 


A{W) 


B{W) 


A{W)B{W) 


2 


2.296 


2.601 


0.853 


2.219 


3 


2.646 


3.096 


0.829 


2.567 


4 


2.792 


3.343 


0.874 


2.922 


5 


2.867 


3.492 


0.834 


3.262 


6 


2.912 


3.586 


0.998 


3.579 


7 


2.941 


3.695 


1.055 


3.898 


8 


2.961 


3.754 


1.118 


4.197 



Table 1-2. From the one-dimensional scaling procedure these values can be extracted. The 
expression of the gap for the infinite long periodic strip is given in equation ( 1.5. 2| ) . 



The combination A{W)B{W) also shows a trend to the two-dimensional situation. 
In two dimensions the critical exponent u < 1. The first derivative with respect to 
the field H, which is here A{W)B{W) and in equation (JjLTli ) it is u{H - H^f-^, 
should diverge at the phase transition point. So 

lim A{W)BiW) = oo. 

W^oo 

The trend is in the data but W = 8 is far to small for this combination to grow 
excessively. 

1.5.2 2-Dimensional, ^ oo , L = oo. 

Standard scaling methods as described in can be implemented now; 



A{h,W) = —g{hW'/',u,Wy^). 



;i.5.3) 



The reduced field h is given hj h = H — (remember He = Hc{oo) ); true two- 
dimensional scaling behaviour is considered, in contrast with equation ( |1.5.1D . The 
finite-size corrections are anticipated by the inclusion of the irrelevant field Ui and 
exponent i/i. In the previous subsection also an expression for the gap A was given 
in equation (|1.5.2|) . As both, this one and the scaling relation above must hold, they 
have to be identical: 

A{W)B{W) [h + H,- H,{W)] = ^g{hW'/\ u,Wy^). 

The left side of this identity is linear in h. The right part also contains higher order 
terms, which we will neglect. Expanding both sides up to linear order yields: 



A{W)B{W)W[H,- H,{W)] 



(7(0,0) + iyJ^»n,|^(0,0), 

OUi 



Aiw)Biw)w = iy'/"^(o,o) + iyi/"+^'Ui-^(o,o). 

oh OUi oh 
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Literature 


first L ^ oo 
then W oo 


L = QW oo 


He 


3.0444 


3.0449 


3.0439 


V 


0.63029 


0.61 


0.62 


Vi 


-0.83 


-1.21 


-1.2lW 



Table 1-3. Comparison of the results for the critical properties of this work with those by 
Blote . The first result corresponds to the two step approach: the length is scaled to 
infinity, L ^ oo, afterwards the width IF — > oo. The second approach is scaling with fixed 
aspect ratio. the value for the irrelevant exponent was taken from the anisotropic scaling. 



The derivatives denote differentiations with respect to the first or the second argu- 
ment. Given A{W), B{W) and Hc{W) for widths W = 2, . . . , 8, both the critical 
properties He and u and the fitting parameters g{0,0), Uidg{0,0)/dui, dg{0,0)/dh 
and Uid^g{0,0)/ duidh can be obtained. 

The results are listed in table The critical field we obtain is of similar quality 
as obtain by cluster Monte Carlo calculation by Blote 0. The critical exponent is 
of substantial less quality. Still, it remains striking to observe that this quality can 
already be achieved by considering such a narrow systems, W < 8 \ 



1.5.3 2-Dimensional, L,W ^ oo with fixed aspect ratio L/W. 

The aspect ratio Q = L/W is fixed and the second step ( |1.5.3|) in the previous 
approach can be implemented immediately; 

A{h,W,L) = jf{hL'/'',u,Ly',Q). 

As we consider fields close to the critical field, \h\ <ti 1, and the corrections to scaling 
are expected to be relatively small, this is expanded up to linear order: 

LA = /(0,0,g) + L^'M,|^(0,0,Q) 

OUi 

+ftLV.(|,0,0,O) + L»...^(0,0,Q)). 

(1.5.4) 

The critical properties H^, v and yi are not dependent on the aspect ratio Q. All oth- 
ers, /(0,0,Q), Mi9/(0,0,Q)/9Mi, 9/(0,0,Q)/9/i and u^d'^ f{^,^,Q) jduidh, do. The 
parameters were extracted by a least square fit to all of the data available. The systems 
considered had a range of properties; widths 4 < < 8, Q = 2,2.5,3,3.5,4,4.5,5 
(the fractional aspect ratio's were only considered with even width W) and fields 
3.00 < H < 3.10. In total 471 points were fitted to this behaviour. An optimal fit 
can be made for a large range of parameters. We therefore had to fix Ui = —1.21; 
the value found in the anisotropic scaling. If the least square fit is made unrestricted, 
unreasonable critical exponents will result. The outcome is listed in table |l-3. 
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1.6 Conclusion 

In this chapter we derived the critical field H^. and exponent v of the two-dimensional 
ITF. Two routes were taken: one was to scale the length of the system to infinity and 
afterwards the width. The other route followed the usual path of scaling; the linear 
dimension was rescaling in both directions. 

The two-step scaling approach was prompted by the strengths and weaknesses of 
the DMRG. Table |1-3| clearly indicates that this makes sense; the accuracy of the 
critical field He may be considered unexpected in view of narrow systems considered. 
The critical exponents are not too impressive. 

Although the essential ingredient in these calculations is the DMRG, no time was 
spent on the procedure or the properties. The next chapter contains that half of the 
research. The complication that was mentioned, strong limitations on the width of 
the system, is a very general feature of the method. Later on, in chapter ^ where the 
frustrated Heisenberg is investigated, the same problem will reappear. 

Apart from two constants. He and z/, this chapter has given us an excellent play- 
ground to test the DMRG and a first taste of finite-size scaling in quantum magnetism. 
This will be of use for the frustrated Heisenberg model. 



2 Critical properties of the ITF through 
DMRG calculations 



2.1 Introduction in the Density Matrix Renormal- 
isation Group 

In this chapter the density matrix renormahsation group (DMRG) is introduced and 
studied. This method was originally proposed by White |^ in 1992 to resolve some 
of the problems that the real space renormahsation group (RSRG) suffers from; in 
contrast with its successful treatment of the Kondo problem, the RSRG gives — 
relatively — poor results for strongly interacting quantum chains. 

It is well established by now, that the DMRG can handle this class of lattice 
problems extremely well, making it the method of choice for one- dimensional quantum 



lattice systems. Already in the initial papers ||5T|] the promise became apparent. The 
treatment of the spin-1 Heisenberg chain demonstrated the capabilities of the 
method. The accuracy was unprecedented. 

The method is variational and systematic. If the computation is scaled up suffi- 
ciently, the approximation the ground state that is made, has to be accurate. The 
surprise lies in the fact that this situation can easily be achieved with a minor com- 
putational effort. Some theoretical investigation on the grounds for this tremendous 
accuracy has been made |^, but no clear understanding exists at the present. At 
the same time the applications have started to diversify. Three main trends can be 
observed. 

The first one is usage of the method in two-dimensional classical problems. Bursill 
et al. followed the standard route of a transfer matrix description to transform a 
two-dimensional classical model into a one- dimensional quantum system. Carlon and 



Drzewihski |T^ build on this extension to settle some of the outstanding issues in 
that field. 

The second diversification is in the direction of chemical molecules. Historically, 
physicists and chemists alike have tried to develop simple models of chemical com- 
pounds that allowed a theoretical treatment. A good example is the Su-Schrieffer- 
Heeger Hamiltonian to model the valance electrons in a polymer. It simplifies the 
individual atoms to lattice points, whereby a one-dimensional quantum lattice model 
is formed. For a proper treatment the Coulomb interaction has to remain long ranged 
and should not be restricted to individual sites. Despite this long-range interaction, 
it seems that the DMRG can still achieve good accuracy although the issue is not 
completely decided yet 0. 

The other direction, upgrading the method to deal with more realistic molecular 



models, has also been tried. White and Martin [56| have applied the DMRG to the 
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orbital description of water. They concluded that the accuracy compares favourably 
with many other numerical methods. 

The third extension is also our main interest, the application to wider systems. 
Noack et al. @] focused on the two-leg Hubbard ladder, while White and Scalapino 



||55[] investigate the t-J model on a 16 x 8 lattice. We were specifically interested 
in the limitations of the DMRG on wide systems. It seems that DMRG is not well 
suited to handle two-dimensional quantum systems, although with tricks and brute- 
force computer power reasonable results can be achieved. In this chapter we explain 
what the difficulties are. 

As a testing ground we use the ITF. The ITF constitutes of a well understood 
and simple case of a two dimensional model with a quantum phase transition. The 
performance of the DMRG is expected to decrease in the vicinity of a phase transition, 
making the ITF our model of choice. Moreover cluster Monte Carlo calculations Q 
have resulted in accurate numerical values for its critical properties to which we can 
relate. 

The aim of the method is to find an approximation to the ground state wave 
function. This is done by 'bootstrapping'; the approximation to the ground state is 
improved iteratively. To represent the ground state a basis in the Hilbert space is 
necessary. This basis is systematically improved by the DMRG. Typically, the route 
taken is to start with a basis for one site and iterative enlarge the Hilbert space by 
adding sites. Without further ado the size of the basis would then grow exponential 
like the Hilbert space at each addition. To avoid this, the most relevant basis states 
are selected and all others are removed. 

The route we follow in this chapter is first to introduce and discuss the density 
matrix which is the key ingredient of the method. Next the geometrical properties 
of the systems we study are listed. Of the two procedures we implemented, one has 
extensively been described in the literature |51|, |19[. The other one was introduced 
by us [|l5l and will be reviewed here. With the method in place, a link with the 
RSRG is made. After that, the performance and the fiexibility of the two procedures 
is compared. The original proposal by White is the more fiexible of the two whereas 
our implementation has the potential to be the fastest. 

The remainder contains two less related sections, one on the actual results for the 
ITF model and one on implementation issues. 



2.2 The Density Matrix 

Suppose we have a state |0o)- This can be the ground state, but for the moment we 
leave that open. In the DMRG we want to find a basis to represent this |0o) "well as 
possible. For the entire system, this statement is trivial; the only required basis vector 
is the state itself. If on the other hand the system is divided in a part A and B, both 



containing a number of sites, as depicted in figure |2-1| , we can find the most relevant 
fraction of the basis in part A. Our aim is to truncate the basis to this relevant part 
and remove the rest of the basis. A restriction to these basis states should allow an 
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A 



B 



Figure 2-1. The system is divided in a left part A and a right part B, both containing a 
number of sites. With both these parts an incomplete basis is associated. In the future, a 
generalisation will often be used; the system is split in three or four parts by treating either 
one or both of the intermediate sites (last site of A and first site of B) separately. 



optimal representation |<^o) of the ground state |(/)o). By optimal we mean that the 
truncation error P, 

2 

I0o) - |0o) , 

is minimal. We have to select those m basis states in A that minimise the truncation 
error P. 

Now that the aim is conceptually clear, we will perform the algebra and obtain 
the optimal set. We have a basis {\i)} in A and {\j)} in B to represent |0o)- Note 
that these bases are not necessarily complete. The only requirement is that the space 
contains |0o), 

I0o) = Y.^ii\^)\j)- 

We want to select the m most relevant basis states from part A. They are 

contained in the basis and orthonormal, thus 

l«) = E<K) ' («!"') = E<*<' = (2-2.1) 

i i 

The best approximation |0o) can be expressed in these states, 

aj 

and the truncation error P, which we want to minimise, is a function of all these 
parameters, 

2 



|0O) - |0O) 



ij \ a / 



(2.2.2) 



The truncation error has to be minimised over the parameters, whilst the orthonor- 
mality condition ( p.2.1| ) is fulfilled. Lagrange multipliers introduced for this 

purpose. The optimum satisfies 

dP dP dP 



duf d4>aj dfi 



0. 
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We introduce the density matrix pai of A, 
pit' = J2^iA'j- 



After some elementary algebra the prefactors (paj and the Lagrange multipliers /i^ 
can be obtained, 



The remaining equations read 

Y^Pii'K - H ufufp,'i»u'^„ = 0. 
i' i'i"a' 

This equation becomes more transparent when we switch to vector notation, 



■ p ■ u" = 0. 



This equation can be interpreted as a projector acting on the vector p- u". This 
projector removes all components of the vector along one of the vectors as 



= 0. 



The basis spanned by {u"} has to contain all vectors |p ■ u"|. A set of eigenvectors 

of p clearly satisfies this condition. Set {u"} to be eigenvectors of p with eigenvalues 
Aq so p • u" = AqU". We know that 



> 0, 



The complete set of eigenvectors u" is orthonormal, so 



II' J a 



ii'j 



(2.2.3) 



To decide which eigenvectors to select, we note that if m vectors are selected with 
eigenvalues A^, then 



P = 1 - ^ A„. 



(2.2.4) 



a=l 



This relation can be derived by inserting and (paj in the definition of the truncation 
error (^^). 
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Given ( p.2.4| ) it is immediately evident that the eigenvectors u" of p corresponding 
to the m largest eigenvalues Aq, have to be selected to build the basis {|q;)}JJLj^. 

It is worth to continue these algebraic manipulations a bit further and show how 
the optimal basis in part A is related to the optimal basis in part B. This can 
be derived step by step. For instance White came to similar conclusions when 
investigating the consequences of a singular value decomposition. Here we simply 
present the result and discuss its consequences. 

The relevant bases states {|q;)J^=i} in B are given by 



) = Y.^'!\j), 



where the vector v" is related to the vector through the relations 

1 



These bases states are orthonormal and they are indeed the largest eigenvectors of 
the density matrix ^jji = J2i (Pij4>ij' as can be seen by insertion; 

1 * / — 

j' V^a iiiji i 

The ground state obtains a very elegant form in these bases, 



l0o) = E'^'.K)b) = E (2-2.5) 

ij a=l 

To every basis state |a) in part A there corresponds exactly one \a) in part B. The 
set size fh is the minimal number of states available in either part A or B. This 
equation (|2.2.5| ) makes ( |2.2.4| ) trivial; the approximation |0o) is simply be given by 



Let us elaborate on the equation ( p. 2. 5 ), as there are a few important consequences 



of it. If there are as many states in B as we wish to select in A, two further assessments 
can be made. 

First, the truncation error vanishes, P = 0; there is no approximation made in 
transforming |0o) to |0o), |0o) = |0o)- AH properties of the wave function in the 
truncated basis are identical to those in the full basis. 

Second, if to every element of the basis in S a set of quantum numbers can be 
assigned, then the distribution of the quantum numbers in part A, contained in the 
set is fixed. Every element has the correct quantum numbers to pair up with 

the quantum numbers of his partner in B to form the required quantum numbers of 
the wave function |0o)- 

In the past few paragraphs we developed a method to distinguish important states 
in a part of the system by the density matrix. This selection criterion can be gener- 
alised to incorporate several wave functions |0/3). This is not a hypothetical situation 
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since we will actually target several states in chapter ^ . A weight Wj3 is introduced 
to differentiate these states {(pp) in importance. We now minimise 

The underlying framework of minimising P remains linear algebra and therefore the 
density matrix becomes the linear superposition of density matrices for each (3\ 

Pii' = J2^/3Pii'- 

An approximation scheme based on the density matrix would be useless if this trun- 
cated basis does not represent |0o) properly. Clearly the indicator to study 



is the truncation error P. Figure (p-2|) shows that the density matrix can indeed be 
an excellent selection criterion as the truncation error P falls off exponentially with 
the number of states kept m. In this case only m = 13 states need to be kept of the 
2^ = 512 to obtain an accuracy of P ~ 10~^ in the wave function. 

In the following procedure we will make approximations to the ground state and 
select the optimal bases to represent them. As these approximations are not identical 
to the ground state, some caution is needed in using the truncation error P as a 
measure of the accuracy; As stated before, the truncation error will vanish, P = 0, if 
the environment contains as many states as we want to select. This does not mean that 
is a perfect representation to the ground state. Instead it means that |0o) = |0o)- 



2.3 Geometry and symmetries of the ITF 

As a last stop before commencing the description of procedures to implement the 
density matrix principle, we have to define the system geometry and investigate the 



consequences on the symmetries of the ITF. As is depicted in figure 1^1 consider 
systems of sizes L x W. The length L is in general larger than the width W. The 
system is periodic in the width direction and open in the length direction. It will 
be split in a left-hand part A and a right-hand part C, both containing m states. A 
intermediate band B, containing the complete basis of 2^ states, separates them. 

The Hamiltonian of such a system contains many symmetries that we can incorpo- 
rate in our calculation. The general form of the included symmetry operators is that 
they are the direct product of three components. Each component acts on one part 
of the system only. For example, consider the translation operator T in the width- 
direction. This operator is the direct product of three translations in the individual 
parts; T = Tj^TbTc- The same holds for the reflection TZ in the same direction,??. = 
TZaR-bTZc, and the spin-reversal operator S = exp^iir ^ij + LW/2) = SaSbSc- 

The ground state \tpo) of the system is translational, reflection and spin- reversal 
invariant; T|?/^o) = T^li^o) = <S\iIjo) = \ipo)- For systems of infinite size in the classical 
ordered region {L,W ^ oo and if <^ 1), it will become degenerate with a state that 
is spin-reversal anti-symmetric. In order to take advantage of the symmetries, the 
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Figure 2-2. The truncation error P as function of the number of states kept m. The system 
is the ITF with periodic boundaries in both directions and L x W = 6 x 3. Part A contains 
the left 3 X 3 = 9 sites. H = 3.0. 
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L 

Figure 2-3. The systems we consider are of dimensions L x W. The system contains three 
parts: a left-hand part A (shaded), an intermediate band B (black) and a right-hand part 
C (white). At every DMRG step part A and B are contracted. 
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bases of part A, B and C are chosen to be eigenvectors of the symmetry operators 
T and S. IZ is used later on. So if {|a)}, {|&)}, {|c)} are the bases of the individual 
parts then 

r^|a) = e^'^^la) , SaW) = Sa\a). (2.3.1) 
Similar relations hold for the other two sets. Thus 
I0o) = 5l0afec|a)|&)|c) 

abc 

and application of the symmetry operations together with (|2.3.1 ) yields: 

ka + kb + kc = Q mod 27r , SaSbSc = 1- 

It is also possible to set up the program to find the lowest state in other symmetry 
classes by forcing other values than and 1 in the equations above. 

The Hamiltonian can be written as the sum of Hamiltonians within the separate 
parts: T-Ca,'Hb and He combined with interactions between parts: 'Hab.'Hbc and 
"^ca; ^ = 'Ha + li.B + 'He + 'Has + 'Hbc + 'He a- To show how to implement the 
symmetries, we will discuss one element of both types. 

First T-Ca'- it is translational and spin- reversal invariant, thus 

{a'lHAla) = {a'\rx^nATA\a) = e''^^^-''-'\a'\HA\a) 
= {a'\Sj^HASA\a) = Sa'Sa{a'\HA\a) 

These relations lead to the appearance of delta-functions, 

{a'\nA\a) = (a'|^yi|a)5s,,,s„4,,,fca- 

It only contains elements within symmetry classes, as one would expect. 

Second Hab- once again, it is translational and spin-reversal invariant. Moreover 
it can be written as 

w 

'HaB = -4 X! ^tn^l+l,n 
n=l 

W 

= -AY.{TArB)-^'Sl,St^,^,{TATBr-^ (2.3.2) 

n=l 

where / is the length of part A. S^j fiips a spin, so S^ jS + SS^j = 0. Inserting this 
and ( gXID in (|2X^ ) gives 



{a'\{b'\'HAB\a)\b) = -AW{a'\Sl,\a){b'\S{'_^,Jb)- 

This substantially reduces the computational effort. Finally: the reflection operator 
TZ is used to make matrix elements like {a'\Sfi\a) real. In fact we could have used this 
last symmetry TZ more, but it only reduces the effort by a factor of 4 while making 
the program far more complex. 
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2.4 Procedures 

In the section on the density matrix, a selection criterion was developed. Now this 
criterion will be incorporated in a procedure. We have employed two distinct proce- 
dures: the first one was originally proposed by White and trivially extended to 



two dimensions by Liang and Pang ^8\. In this scheme, part A, figure is enlarged 



iteratively at the expense of part B. Each step one site is moved from part B to part 
A. The basis on part A will thus incorporate more and more sites, while the total 
system size remains the same. Excellent introductions already exist to which we refer 



5T| , |T9[. The second procedure we introduced ourselves in This will be reviewed 
next. 

When the split-up of figure is made, it is tempting to use the ID DMRG 
method directly: a site is replaced by a band. The ground state |0o) of the entire 
system ABC is calculated and the optimal basis for block AB is selected through the 
density matrix. However, one runs in severe difficulties. In the section on the density 
matrix it was shown that when we select as many states from part A combined with 
B as there are in part C, the outcoming basis has to have the appropriate quantum 
numbers be to combined with part C. So if we also would transform part C in the 
density matrix basis, we can write 

I0o) = J2 \/>^\a)\a). 

a 

We know that T|0o) = S\(f)o) = \(po), thus 

Tc\a) = e*^"|a) , Sc\a) = Sa\a). 
This will immediately dictate the quantum number of the newly built 

7^ Tela) = e~^^"\a) , SaSbI'^) = Sala)- 

Thus Sa = Sa and ka = —ka- The distribution over the symmetry classes in part C 
forces the selected states in block AB to be in "conjugate" -classes. To overcome this 
problem, we need to increase the number of states in part C. In that case we can 
really make a selection and it allows us to change between symmetry classes. 

In the ID procedure the solution is to add one extra site to the environment. 
The number of states in the environment is then doubled. In our set-up this would 
correspond to adding an extra band between B and C. This is computational far too 
expensive. 



We now introduce variants on White's infinite-size and finite-size algorithms ^1 
that increase the number of states in the part C. 

First we consider our infinite size approach. We only have to describe one step in 
the process as it is an inductive method. We have a basis of m states for a system of 
length I. We then proceed as follows: 



We construct the combined system as depicted in figure p-4| -a by taking this 
basis in part A and C together with the complete basis in the intermediate 
band B. {L = 21 + 1) 
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We calculate the ground state \(f)Q) and obtain m basis states for a system of 
length / + 1 by orthonormalising {|/3c)}, 

l/^c) =J2^o,bc\a)\b), 



and span the space we previously denoted by {| 

Suppose that block AB has / symmetry classes. To every symmetry class we 
add m/ f basis states constructed randomly from the m2^ states in A and B. 
We end up with m + f ■ m/ f = 2m basis states for a system of length I + 1. 

In part A we now take the m basis states for a system of length / and in part 
C we take the newly constructed 2m states for length / + 1. (L = 2/ + 2) This 
yields the configuration in figure ^-4|-b. 



• We calculate the ground state |0o) and obtain 2m basis states for length / + 1 
by orthonormalising {|/3c)}- We replace the basis of part C by this basis and 
repeat this step a couple of times (~ 3). 

• We select from the 2m basis states for length I + 1 m states on basis of the 
density matrix. 

Now we have returned to the original situation with the exception that / has increased 
by one. The new ingredient is thus to add m random states to the basis and iterate 
until the result has converged. 

In the same line our finite size approach lies. Suppose we have basis sets of m 
states for lengths l,L — I — 1 and L — I — 2, where L is now fixed and independent of 
/. The iteration step consists of the following actions: 

• We take the basis for / in part A, the basis for L — / — 1 in part C and the 



complete basis of the band in part B. See figure 2-5 -a. 



We calculate the ground state |0o) and obtain a basis for length Z + 1 by or- 
thonormalising {|/3c)}- 

In the same way as in the infinite-size algorithm we add m randomly chosen 
states to this basis. 

In part C we take the 2m basis states for length / + 1 and in part A the m 
states for length L — I — 2. This is depicted in the first of the two pictures in 



figure 2-5 -b. 



We calculate the ground state |0o) and obtain 2m basis states for L — I — 1. 



In part C we take the 2m basis states for length L — / — 1 and in part A the m 
states for length Z; see the second picture in figure |2-5[ -b. 
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Figure 2-4. A inductive step in the infinite-size procedure consists of a start-up to obtain 
an initial approximation for states in a system of length I + 1 and iterative calculations 
to make the basis converge. The numbers in the rectangle are the number of states in the 
parts. The intermediate band B always contains the complete basis of 2^ states. 
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Figure 2-5. A inductive step in the finite-size procedure also consists of a start-up to obtain 
an initial approximation for states in a system of length / + 1. Afterwards we move back 
and forth between lengths / and / -|- 1 to make this converge. 



• We calculate the ground state |0o) and obtain 2m basis states for I + 1. These 
last four steps are repeated a couple of times (~ 3) 

• We select from the 2m basis states for length I + 1 m states on basis of the 
density matrix. 

Once again we have returned to our starting position while increasing the length / 
by one. By sweeping through the system we can therefore systematically improve 
the basis. This method convergences at a similar speed as the ID approach; after 3 
sweeps through the system the final result is achieved. 

The approximation scheme to the ground state is both variational and systematic. 
The energy of each state in the selected subspace is higher than the one of the ground 
state. As this includes the state we select by minimisation of the energy, our estimate 
to the ground state energy is variational. In section ( [4.2| ) we will meet other examples 
of variational principles that can be used within the DMRG scheme. 



2.5 Connection with the Renormahsation Group 
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The systematics come in from the iterations. When a state is selected, it is trun- 
cated and transformed along with the basis. The next iteration it is used as a starting 
point for the minimisation routine. If the density matrix eigenvalues drop off fast 
enough and the truncation thus does not severely alter the state, the energy of the 
starting point will be almost the same as the outcome of the last iteration and the 
during the minimisation this estimate can improve further. 

After several sweeps through the system, the energy will start to oscillate. This 
effect is small, but the cause of it gives ground for further improvement. Suppose we 
have found a state |0o) at the previous iteration. This state is truncation to |0o). 
In the truncation a part of the state is lost, so |0o) 7^ |0o)- The truncated space of 
|0o) lies in the subspace of |0o) though. We know that |0o) corresponds to the lowest 
energy in that subspace, thus 

(0O|H|0O) ^ (0o|^|0o) 



(0o|0o) (0o|0o) 

|0o) is transformed and thereby embedded in an other subspace. It is used as a starting 
point for the next minimisation. The next estimate |0o) will thus have a lower energy, 

(0O|^|0O) ^ (0O|^|0O) 
(0o|0o) ~ (0o|0o) 

No strict statement can be made on the relation between the energies of |0o) and 
|0o)- The energy does not need to decrease monotonically and will start to oscillate. 

Once this happens, we can assume that the basis states we use, lie in the most 
relevant symmetry classes and it is no longer necessary to increase at every location 
in the system the number of states in part C from m to 2m. Leaving out the iterations 
to increase the number of states in part C, the previous estimate |0o) can still be 
represented exactly after the basis truncation , |0o) = |0o)- This is explained in the 



paragraph after equation (|2.2.5| ). The ground state |0o) of the next step satisfies 



(0o|7^|0o) _ (0o|7^|0o) ^ (0o|7^|0o) 

(0o|0o) (0o|0o) ~ (0o|0o) 

The energy will now decrease monotonically. 



2.5 Connection with the Renormahsation Group 

The name DMRG has led many people to believe that it is another implementation of 
Wilson's renormahsation group. This is incorrect in at least two respects. First, this 
method does not contain a semigroup operation; where in the renormahsation group 
the Hamiltonian is mapped onto itself with different parameters, here we change — 
read truncate — the Hamiltonian all together. The name group is unfortunate and 
misleading, but it is well established and to avoid further confusion, we will stick to 
the name DMRG. 
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Figure 2-6. On the left the wave function of a particle in a box is depicted. At current a 
continuous version of the DMRG does not exist, so one should consider the location in the 
box a very fine grid of discrete points. On the right we see the selections of the DMRG, 
one state suffices, and of the ordinary renormalisation group (RG) where many states are 
needed. 

The second distinction lies in the fundamental difference of basis selection. The 
DMRG first calculates the approximate ground state of the entire system and af- 
terwards makes a selection via the density matrix of this state. The renormalisation 
group, on the other hand, does not consider the properties of the entire system. The 
basis states are solely selected on merit of the energy within their own part of the 
system. 

The well-known example [0, Q of a particle in a box with impenetrable walls 
illustrates this selection criterion clearly. By its graphical simplicity figure p-6| demon- 
strates that in the DMRG one basis state suffices whereas the renormalisation group 
needs many to represent the wave function properly. 

2.6 Performance 

In this section a discussion is given on the accuracy that can be achieved by the 
DMRG procedures. First we consider issues that affect both methods, afterwards a 
comparison between the two methods is made. 

2.6.1 General limitations. 

The most striking limitation is that only systems of small width can be handled. 
As the width W increases, the number of states m kept in the procedure has to 
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Figure 2-7. The accuracy of the DMRG method for given number of states m (numbers in 
graph) as function of the width W. H = 3 and L = 20. The system is periodic in both 
directions. The reference value is taken from a DMRG calculation with m = 128. 



increase exponentially to maintain the accuracy |28|. In figure |2-7| this is exemplified 



for the ITF. To get a flavour of the background to this behaviour, we will discuss a 
pathological example where this statement can be proven exactly. 

Consider a system containing a set of chains as depicted in figure P-8| . For each 
chain a DMRG calculation can be performed and it is found that rriQ states are 
required to obtain a given accuracy. To achieve the same accuracy for this system of 
width W = 4, for each chain tuq states have to be preserved as there is no interaction 
between the chains. The basis of the entire part A is a product of bases for the 
individual chain pieces. Therefore there are m = states necessary for the same 
accuracy; a clear proof of the exponential growth of the number of states m with 
increasing width W. 

In the study of the ITF we are fortunate enough to go somewhat further than 
this rough-and-ready argument. Far from the phase transition {H ^ 1 or ^ 1) a 
connection can be made with perturbation theory. 

Consider the quantum disordered phase on a periodic system. Split the Hamilto- 
nian into a unperturbed part Ho = 2if and a perturbation V = — 4 X]j j Sf^^Sfj^^ 

S^j^i). We split the periodic, rectangular system of size L x W again in two parts; 
A and B of sizes I x W and {L — I) x W where / is an arbitrary length smaller 
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Figure 2-8. The system contains only interaction within the chains. Between the chains 
no interaction exists. The system is spht up, and the DMRG procedures grow part A on 
expense of part B. Either site by site or row by row. 



than L. They both contain 2W spins that border the other part. The unperturbed 
ground state |0) has all spins pointing down in the «S^-direction. It is the direct prod- 
uct of two equivalent states restricted to A and B; |0) = |0)yi|0)s. We know that 
7^0 1 0) = -HLW\0) ^ Eo\Q). Perturbation theory yields 



(2.6.1) 



The perturbation flips a pair of neighbouring spins. This pair can be in a single part 
or it can cross the border between both parts. In the latter case the spins are adjacent 
across the boundary between part A and B. Deflne {la)^} to be the set of states with 
the flipped pair in part A. Analogous for Moreover let be the set with 

one spin flipped on the nth boundary site with B and define in an equivalent manner 
The perturbation expansion can now be rewritten 



= |0)a|0)s + ^ (e I«)a|0)b + E I0)^l^)^ + E \n)A\n)B) 

\ a b Ti / 



+0 



(2.6.2) 
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As if 3> 1, it is necessary to reproduce all these terms for an accuracy which is 
equivalent to the first order perturbation theory. The minimal number of states needed 
in part A is therefore 1 for the first term in ( ^.6.21 ) plus 2W for all the boundary terms. 
We have confirmed this prediction explicitly in both the small and large H limit ( 
H = 1/50,50). 

The same line of reasoning also holds for the second and higher order perturbation 
terms. We expect for an error comparable to the nth order perturbation theory that 
m ~ W"', SE ~ This is always an upper bound for number of states m 

needed, m < W"' for a given accuracy 6E ~ Only when the different orders 

in perturbation theory become distinguishable in size — the limit of large H — the 
equivalence holds. Through combinatorics even the prefactors can be calculated. 

Both arguments above indicate that it is of the utmost importance to limit the 
interaction between the parts, although the statements are not identical. If there 
exists a representation in which the interaction between the parts is still large, but 
where the perturbation expansion needs only few states, the DMRG easily obtains 
highly accurate results within that representation. 

The question what basis to start from is still open. Xiang chooses to treat 
the Hubbard model in momentum space. This representation introduces an extra 
quantum number that can be conserved: the momentum. As this leads to a relatively 
strong restriction on those basis states of the different parts that can combine, the 
number of states kept m could be increased substantial. Moreover we expect to have 
a better starting point for the ground state by this conservation. On the other hand, 
from the view of perturbation theory it is unclear whether the orders in perturbation 
theory can be reproduced easily as the interaction between the parts becomes much 
more complex and much more extended. 

In the case of fermions it is easy to change the representation of the particle 
creation and annihilation operators from real space to for instance momentum space. 
It is also possible to switch to the mean field quasi particle representation in which 
the ground state contains no quasi particles. This is done in the BCS theory by means 
of a Bogoliubov transformation. For a DMRG calculation such a representation may 
also be not ideal despite the fact that the mean field ground state can be represented 
by exactly one basis vector (with no quasi particles present). It could be a source of 
bias in the final result. 

Other attempts have been made to avoid the issue of limited width all together by 



smoothing the boundary conditions and thereby reducing the boundary effects p9 



Another restriction is related to phase transitions. A clear indication of the vicinity 
of a phase transition is divergence of the correlation length. At the phase transition 
algebraic behaviour of the correlation functions is expected. This is conflicting with 
the manner in which the ground state is built by DMRG. Successive basis rotations 
and truncations where at every step a new site is included in the basis clearly favour 
exponentially decaying correlation functions ||35|. For the ITF the accuracy of the 
DMRG at the phase transition indeed decreases although it remains unclear whether 
this is a general feature of phase transitions. Figure |2-9| illustrates this. 
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Figure 2-9. The accuracy of the DMRG method for different number of states m (numbers in 
graph) as function of the field H. The system is periodic in both directions with dimensions 
W = 4: and L = 20. The reference value is taken from a DMRG calculation with m = 64. 
The critical field of the two-dimensional system is He = 3.044 

Several basis rotations and truncations applied to operators cause them to deteri- 
orate. For the intermediate band or site the basis is complete and the representation 
of the operators on these sites is perfect, but for sites further away from the part 
boundaries, some basis truncation have been applied and the representation of the 
operators has lost the connection with the states truncated from the basis. A con- 
sequence of this is that the correlation functions for sites that are simultaneously 
or successively included in a part by the procedure are of higher quality than those 
between sites that lie further apart. 

2.6.2 Comparison of both methods. 

In section |2.4| a version of the DMRG was introduced that added bands to a part. 
White's original proposal was to add sites to a part. These two procedures, site and 
band, can be compared on two grounds: 

• Computational effort with given accuracy. 

• Flexibility of the procedure. 
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Let us address these criteria in the same order. 

A straightforward comparison of the two computer programs favours the site 
method heavily. Unfortunately this says more about the software than about the 
quality of the method. A fair approach contains two steps: first the accuracy as 
function of the number of states kept m is related. Afterwards the computational 
effort of both procedures as function of the number of states m is estimated. We have 
done both for the ITF. 

Figure |2-10| depicts the first step. The band method clearly needs fewer states for 
the same accuracy. 

The bottle neck in the computational effort is finding the ground state |0o) in a 
given subspace. This is done by applying the Hamiltonian in the order of 5 times to 
the wave vector. In the site procedure this costs ■ operators per projection. 

The factor 1/2^ follows from conservation of the spin reversal operator S. To sweep 
through the system there are thus tgue = ^-m? ■L-W'^ operations needed. The band 
method has a far larger space in part B but it is possible to use translational invariance 
in the width direction and fewer steps have to be taken to sweep through the system. 
Per step it costs m? ■2'^ /W operations (m > 2^). With the moving back and forward 
at one location (3 successive diagonalisations) , we get t^and = 3 ■ 5 ■ rn? ■ L ■ 2^ /W 
operations for one complete sweep. Up to width W ^ 8 these two estimates are 
similar, 

tsite 5m^LWy2'^ 1 



hand lbm^L2^/W 3 2^^+2' 

with the band method achieving higher accuracy with the same number of states m. 
After that the site method clearly becomes faster. 

To summarise the comparison of the computational performance: as long as the 
systems do not become too wide {W < 8 ) , both methods are comparable in efficiency. 
The band procedure is slightly preferential as it could be faster and the ground state 
is more symmetric. The difference in calculation time we observed is mainly the result 
of the quality of the software and specifically the use of the BLAS 0] routines. The 
difference in calculative performance is not conclusive. 

Next we compare the flexibility. The computational performance is in practical 
situations irrelevant as a large fraction of the effort spent in research is dedicated 
to probing properties of models by changing parameters and such. This requires a 
degree of flexibility in the procedure. By this criterion the site procedure is clearly 
preferential as both the geometry and interactions can easily be changed. The band 
procedure obtains its better calculative performance from the translation invariance 
of the ground state. This fixes both the geometry and the interactions to support 
translational invariance. The only advantage of the band method in the terms of 
flexibility is that by strictly conserving the translational symmetry, also states with 
difference momentum can be targeted. 

The overall conclusion is that if one wants to test a model the site procedure is 
preferable, whereas if maximum accuracy is required the band procedure should be 
used. 
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Figure 2-10. The accuracy of the site-procedure (asterisk '*') is compared to that of the 
band-procedure (cross x.) for various system widths W and number of states m (numbers 
in the graph). The transverse field H = 3 and the system length L = 20. The geometry of a 
cylinder is used; Periodic boundary conditions in the short width-direction and open ones 
in the long length-direction. 
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2.7 Some considerations on the implementation 

Choosing the right software language and tools is very much dependent on both the 
hardware and software available. As the progress in these is rather swift, remarks in 
this section will be relatively soon outdated. Still it can be of tremendous use having 
an idea of what the implementation could look like. We state our comments going 
from more conceptual to language specific ones. 

In implementing the Hamiltonian, operators of the form Si ■ S2 have to be ap- 
plied to a wave function It is recommended to apply the members of this pair 
successively, Si ■ (1S2I0)), not simultaneously, {Si ■ S2)\(p)- Algebraically these two 
multiplications are identical, but computational the second one is far more expensive 
than the first. 

Many sophisticated matrix manipulation libraries exist and these can be used. 
For that purpose the intermediate sites have to be merged with the larger parts; in 
the site procedure, the left intermediate site can be merged with the left part and the 
right intermediate site with the right part. In this fashion we obtain two bases, one 
for the left part, {\i)}, and one for the right part, Now a wave function |0) can 

be expressed by a matrix $, 

The use of this representation touched immediately at another implementation 
matter; applying operators to this wave function now transforms into matrix mul- 
tiplications. There are quite a few libraries available that are highly optimised to 
perform this task of matrix multiplication, most notably BLAS 0. 

To end our list of remarks, the programming language that we used is C-I--I-. There 
are two reasons to do this: 

• At current, it is one of the most widely used languages. This ensures the avail- 
ability of tools, libraries and sophisticated compilers. 

• C++ allows object oriented programming. This requires a different way of 
analysing a problem than the well-establish procedural approach, but it makes 
the code more transparent and reusable. Especially in combination with pack- 
ages like STL [| this provides a powerful programming environment. 

Naturally C++ is not unique in either of these properties. FORTRAN90 also fits at 
least the second argument, but is less widely used. 

Object oriented code in C++ is slower than FORTRAN or FORTRAN90, but 
in a DMRG implementation by far the most time is spend on finding the ground 
state. As this is code-wise only a fraction of the entire program, this can easily be 
implemented in a procedural fashion using fast libraries (BLAS). 

^standard template library, part of the C++ library 
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2.8 Results for the ITF 

The scaling analysis of the previous chapter is based on DMRG calculations. The 
ground state of the ITF model lies in the space with spin reversal quantum number 
S = 1. The first excitation appears when one spin is flipped, so S = —1 for that 
state. We calculate in both spaces the ground state energy Eq^^'^ and the gap A is 
given by 

A = Et-' - Ef. 

The gap is thus the difference between two large numbers. This causes the relative 
error to increase by two orders of magnitude and the highest achievable accuracy is 
needed. 



Experimentally it has been observed [54] that the truncation error P becomes 



proportional to the error in the energy. This is not surprising as the truncation error 
P is a direct indicator of quality of the basis truncation. The estimate of the energy 
is improved by extrapolating to zero truncation error, P = 0. 

As we mentioned earlier in section ( |2.6.1| ), the accuracy is highly dependent on 



the length of the boundary between the different parts. In a periodic system, each 
part is coupled on both sides to the rest of the system, so then the boundary is twice 
as long as in a open system. With the accuracy already so much under pressure, we 
cannot afford to double the length of the boundary. 

Unfortunately it is not possible to study system of larger width than W = 8 
with the DMRG. This means that the DMRG method has strong competition from 
cluster Monte Carlo algorithms 0. It has to be stressed though, that in cases where 
the Monte Carlo methods fail as a consequence of the sign-problem, DMRG can 
provide an excellent replacement. 



3 Spin stiffness and finite-size scaling of 
the frustrated Heisenberg Model 



3.1 Introduction 



The interest in the frustrated Heisenberg model arose with the discovery of high 
Tc-super-conductivity. The ground state of the undoped compound shows long-range 
antiferromagnetic order, which can be well described by the Heisenberg model. How- 
ever introduction of a small number of holes destroys the long range order and triggers 
superconductivity. In a way doping and frustration have similar effects on the long 
range magnetic order. The underlying similarity is not really understood but efforts 
are made to connect them, see for instance the book by Auerbach |0. To improve 
the understanding of these effects, it is of interest to study a frustrated Heisenberg 
model. 

In more recent years it has become a research topic of its own. It may serve as an 
example of a system exhibiting a quantum phase transition going from Neel order to 
collinear order with increasing frustration. Hopes are that new and so far unknown 
phases will be discovered in the intermediate regime. 

A link can been made between this frustrated Heisenberg model and the material 
CaV^iOg. If the spins are positioned on a rectangular lattice with every fifth site 
unoccupied, the model may describe CaViOg 



Our intention here is to employ the spin stiffness as a means of measuring the 
order present in the system. In this chapter, the spin stiffness will be calculated in 



the Schwinger-Boson Mean Field (SBMF) theory This mean field theory 

introduces the mean fields in the interactions between spins and allows for the ground 
state to be rotationally invariant. Afterwards the spin length is no longer conserved 



to be ^ and we will discuss the consequences in some detail at the end of section 



The spin stiffness in the SBMF serves as a guideline for the effects of finite- 
size scaling, which is an indispensable tool for the numerical calculations using the 
DMRG method. The next chapters contain the numerical parts of our investigations; 
in chapter ^ the DMRG method is combined with finite-size scaling to obtain the 
spin stiffness of the two-dimensional system. In chapter |^ the correlation function are 
studied by a combination of the DMRG and Green Function Monte Carlo simulations. 
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3.2 The frustrated Heisenberg Model 



This model describes interacting quantum spins on distinct lattice sites. The spins 
Sj have length i. Hence |4f = |. The Hamiltonian is given by 



H = Ji^Si ■ Sj + J2^Si ■ Sj. 

It incorporates interactions between nearest-neighbour pairs (ij) with strength Ji and 
next-nearest-neighbour pairs [ij] with strength J2. We will focus on a square lattice 
in two dimensions with L x W = N sites. At the end of this chapter we will consider 
systems where L ^ W. The next-nearest-neighbouring sites are connected through 
the diagonals of the lattice. The lattice is bipartite; it can be split into equivalent 
sublattices A and B. If site j is in sublattice A then all its nearest-neighbours are in 
sublattice B. The function (—I)-' will therefore be 1 if j G A and — 1 if j G -B. 

Both Ji and J2 are taken to be non- negative. If J2/J1 is small the dominating Ji 
term describes the usual antiferromagnetic interaction and the ground state will be 
Neel ordered. For instance spin-up on sublattice A and spin-down on sublattice B or 
{—iy{Sj) = rus. For the opposite case, J2/J1 large, the system decomposes in two 
Neel ordered sublattices which, however, have the same quantisation axis. This is the 
so-called coUinear ordering. In this ordering alternating strips of up and down-spins 
will occur. Suppose the strips are oriented along the x-axis, then (— = 
where the j-th spin in on the y-th row. 

Clearly these couplings frustrate each other on a square lattice when they are 
comparable in size. It is beforehand unclear what will happen then. Whether a dif- 
ferent phase exists in the region between the Neel ordering and the collinear ordering 
remains uncertain. In the literature there have been speculations ranging from dimer 



- to disordered phases |]16|, ^ ^ In all of these cases it is conjectured that the in- 
termediate phase has no long-range order. The correlation length thus becomes finite 
and a gap in the energy spectrum opens up. This energy gap is a good indicator for 
the existence of an intermediate phase. Unfortunately, it is very hard to achieve high 
enough accuracy in the numerics as it involves the subtraction of two 0{N) numbers 
(the ground state energy is subtracted from the energy of the first excitation) to yield 
a C(l) number, the gap. 

Another indicator for long-range order is the spin stiffness. This quantity poten- 
tially allows for higher numerical accuracy. In the current chapter the spin stiffness 
of both the Neel and the collinear phase will be derived in the Schwinger-boson mean 
field approximation. Part of the material in the chapter is published in in . 



3.3 The spin stiffness 

To investigate the phase diagram, one could look for order parameters. However, if 
the frustration increases, a new and different ordering might appear with an unknown 
order parameter. So, many different order parameters would have to be examined. 
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There exists a far more general approach. It is known that when a continuous symme- 
try is broken, a Goldstone mode appears. This is a 'symmetry-restoring' excitation. 
For this kind of spin model, these are called spin waves, as they are wave-like fluc- 
tuation in the spin orientations. In ||23| it was shown that for both types of order 
considered here (Neel and coUinear) the dispersion relation is linear {u = c|q|) for 
low energies. The velocity c of this mode satisfies 

(3.3.1) 

Here ps stands for the spin stiffness and x± is the magnetic susceptibility of the 
system perpendicular to the orientation of the ordering. A positive spin stiffness 
{ps > 0) is an indication that a broken-symmetry is present, although the actual 
order parameter might be unknown! Therefore it is an excellent measure of magnetic 
order that arises from a broken symmetry. For the Heisenberg model an order related 
to a local organisation of spins in the intermediate range of frustration has been 
suggested |]TB|, The global rotational invariance is not broken, therefore 




the spin stiffness should be zero, ps = 0. In short: the spin stiffness is an excellent 
indicator to distinguish between a broken symmetry ground state {ps 7^ 0) and other 
ground states {ps = 0). 

3.3.1 Direct derivation at zero temperature 



Fisher, Barber and Jasnow |18] showed that the spin stiffness ps, that is related to 
the dynamics of the system, can be obtained by a static twist in the order parameter. 
They imposed boundary conditions to achieve this twist. Instead of following them 



in detail, we will directly focus on the order parameter. In this route has also 
been followed and we could reduce the formula there to zero temperature, but the 
transition from finite temperature to zero temperature is subtle and we will instead 
make an explicit derivation of pg. In this way the concepts involved can be made 
more transparent and subtleties concerning the necessity of periodic boundaries can 
be addressed more directly. The resulting expression for the stiffness appears in several 



places in the literature, e.g. 116 



We consider the case of Neel order as an example. Afterwards it will be shown that 
the expression we found holds in more general cases. Before we start let us mention 
that a two-dimensional system exhibiting long-range order at zero temperature is not 
conflicting with the Mermin- Wagner theorem PO]. The theorem only forbids long- 
range order at finite temperature. 

The ground state of an antiferromagnet satisfies 

(0|7^|0) = Eo, 

(0|5J|0) = {-lyrris. (3.3.2) 

It is invariant under translations over lattice vectors respecting the bipartite breakup 
of the lattice. Later on we will come back to this feature. The search for the spin 
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stiffness ps can be formulated as finding an excitation |q) such that 

(q|7^|q) = Eo+iiVp.|q|' + ..., 
(q|S'J|q) = (-l)^mscos(q- rj), 
(ql^Jlq) = (-l)^m, sin(q- Tj). 

The implementation of a twisted orientation like this is a hard task. It can be avoided 
by a similarity transformation. Define the unitary operator W(q) by 

W(q)=exp|^zq-^r,5|j. (3.3.3) 

It corresponds to a rotation about the iS^-axis. It is easy to derive that 

W(q)5;wt(q) = ^/e^'^-^^ , U{c^)SjU\ci) = Sr^r'^-^^. 

We will now transfer the twist from the order parameter to the Hamiltonian with 
this operator W(q). Define 

H(q) = U{q)nU\q) 

^ [ij] 

and |q) = W(q)|q). Combining the expressions for |q) and H(q) it is trivial to derive 
that this new state |q) satisfies 

(q|^(q)|q) = E, + ^NpM\' + 

(q|5J|q) = (-l)^m,. (3.3.5) 

|q) clearly is an excitation of 7i(q) as the ground state is trivially given by W(q)|0). 
We know that 

lini |q) = |0), 

which can be read as an invitation to apply perturbation theory. Still there is a quite 
subtle issue that must not be overlooked. If we would apply perturbation theory with- 
out further ado, the resulting wave function would be W(q)|0). The order parameter 
would also be twisted and consequently it would not satisfy (|3.3.5|) ; 



{0\U\q)SJU{q)\0) = (-l)^m,cos(q-r,), 
(0|Wt(q)5jW(q)|0) = -(-l)^m,sin(q-r,). 
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The underlying physical idea is that if we apply a twist with wave vector q by slowly 
turning up this q, the system and the order parameter will follow adiabatically and no 
energy increase will occur. This is not what we set out to achieve. We are considering 
the case where an integer twist over the entire system exists. The wave vectors q 
that are allowed then fit the lattice; all terms exp(iq ■ (r^ — rj) can by simplified to 
exp{ic[-rs) where is the smallest connecting lattice vector. Even when and Tj are 
adjacent across the periodic boundary of the system this replacement can be done as 
q ■ (rj — Tj) = q ■ mod 27r. To summarise, we seek the homogeneous state |q) given 
an integer twist over the lattice. Respecting this subtlety, |q) can be derived. 
Define 



J 



q=0 



^ (ij) 

+ ^E(r.-r.)(5;5, 



t = - 



dq 



q=0 



^ {ij) 
J2 \ ^ 

^ [ij] 



A note on the terms rj— is in place here. The factors exp(iq-(ri— rj)) as appearing in 
?i(q) are periodic in the difference rj — rj. Adding or subtracting a vector {nL, mW) 
leaves this term invariant. When taking the derivative it is therefore necessary to 
insert the smallest possible argument for Ti—Vj. So when and Tj are on opposite sides 
of the -periodic- boundary we have to replace them by r^; the connecting elementary 
lattice vector. 

The state |q) is now given in first order perturbation theory by 
1 



q) = + 



7J + --- |0). 



Eq — H" 

This state clearly satisfies the condition on the order parameter (|3.3.5|) as both the 
current- current correlation j and the kinetic term t are translational invariant. The 
expression for the stiffness ps can now readily be obtained; 

I = -4(o|tlo) + ^(o|J-l-/|o) 



= T + J. 



(3.3.6) 



The spin stiffness is a tensor and only when the system has quadratic symmetry 
it can be denoted by a scalar. 

The ground state |0) is reflection symmetric whereas the current j is not. Therefore 
(0|j|0) = and we can safely write the current-current correlation as 



■|0N_ 1 y^ (o|jk)Hj|o) 



(3.3.7) 
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For expression ( |3.3.6| ) to hold, it is essential that the system is translational in- 
variant in the direction of the twist q. Recapturing the subtlety in the derivation, it 
can be readily seen that open boundary condition will not lead to the correct result. 
In that geometry we will find |q) = W(q)|0). We can borrow a physical argument 
from the realm of superfluidity to put this in a broader context. 

Periodic boundary conditions can easily be achieved for a bucket of helium by 
bending this bucket round a massive cylinder and connecting both ends together on 
the other side of the cylinder. The twist we apply corresponds to a steady rotation of 
this cylinder with the bucket around it. The superfiuid will not respond to this twist 
whereas the normal component of the helium will start to rotate along. In the frame 
of the cylinder it is the superfiuid that flows, whereas the normal helium remains 
inert. Consequently, the kinetic energy increases by the superfiuid current. 

Open boundaries correspond to helium in a closed bucket. The walls of a closed 
bucket prevent a current from running and a zero increase in energy will be the result. 
It is worthwhile mentioning that Pollock and Ceperley ||3^ have implemented this 
idea directly in a Green Function Monte Carlo calculation. They follow a particle 
winding around a periodic cell. If it winds faster round a periodic cell than expected 
by Brownian motion, superfluidity exists. 

The formula (|3.3.6| ) is applicable to systems exhibiting all kinds of long-range 
magnetic order. If another type of ordering was twisted with an order parameter 
which was also defined locally, the same route would have led to this general result. 
The assumption of antiferromagnetic order can therefore be dropped and we end up 
with a fairly general expression for at zero temperature. The only restriction is 
naturally that the order must break the symmetry; otherwise Goldstone modes or 
spin waves will not appear and no twist is possible. Equation (|3.3.6| ) is completely 
equivalent with the expression found in [|12| and [llH . 



3.4 Schwinger-Boson Mean Field Approximation 

The behaviour of the spin stiffness in a two-dimensional system can be obtained by 
finite-size scaling analysis of relatively small systems. It is then necessary to establish 
the size dependence of the properties and in the current and following sections we 
will derive these in the SBMF approximation. 

The Schwinger-boson representation is in itself a mere reformulation of the prob- 
lem. The only extra condition that has to be imposed is that the number of bosons 
per site has to be fixed to 1. Any operator correctly transformed into the language 
of Schwinger-bosons will conserve this property. 

Standard mean field theory is based on the assumption that there are only small 
fluctuations in the orientation of the spins. In our case that is a poor approximation; 
quantum fluctuations play a major role. In the unfrustrated case, J2 = 0, Monte Carlo 
calculations [^, ^ have shown that the staggered magnetisation is about 60% of 



the classical value. 

A possible improvement of the mean fleld is to incorporate correlations between 
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neighbouring spins. The Schwinger-boson mean field approximation is a first step in 
this direction, it was originally introduced by Arovas and Auerbach [Q] Here we will 
apply this approximation and derive expressions for the energies and wave functions 
of all states of the frustrated Heisenberg model. This will be done both for Neel and 
collinear ordered ground state. 

The first step in this method is to represent every spin by two bosons; 

S+ = a^b, 
S- = ah\ 

S' = ^{a^a-b^b). (3.4.1) 

We restrict ourselves to the subspace in which a^a + 6^6 = 25*, where S is the total 
length of the spin. In this case S = ^. Two observations can be made: 

• The combinations of bosons in ( p. 4.1 ) satisfy all commutation relations of the 



original operators. 

• None of the operators S~^, S~ and will connect this subspace to subspaces of 
other local spin since they all commute with a^a+b^b. Therefore the Hamiltonian 
will only map states in this subspace onto other states in the same subspace. 

The Schwinger-Boson representation is therefore just a reformulation of the problem. 



The transformation in ( |3.4.1| ) does not favour any specific orientation of the mag- 
netic order and the mean field approximation for these bosons will also preserve this 
symmetry. 

Now the mean field approximation can be developed. Many authors have done 
this before and we follow the route taken by Mila et al. |^ for this particular 
problem. As stated before we will consider two types of ordering: Neel and collinear. 

We begin with the Neel ordering. This is expected to be the preferred ordering 
when J2 ^ Ji- 

A first step towards the mean field approximation is to rotate the spins over tt 
around the 5^-axis on one of the two sublattices; 

St - (-1)^5+, 

s. ^ {-irs;, (3.4.2) 
Sf Sf. 



This transformations is inspired by the Marshall sign rule p9|, which states that in 
this new basis the ground state of the unfrustrated antiferromagnet has only positive 
definite coefficients. It is of course not strictly necessary, but it avoids the need of 
complex numbers which otherwise would appear. 

For the mean fields only those combinations of creation and annihilation operators 
can be taken that are invariant under rotation round the iS^-axis. Define 

Vij = a,a] + hb], (3.4.3) 
Bij = Qibj + bittj. (3.4.4) 
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The Hamihonian becomes 

n ■ 



3. 
2' 



where we have inserted ( |3.4.3| ) for the ferromagnetic oriented pairs of spins and ( |3.4.4| ) 
for the antiferromagnetic ones. This expression is related to the fact that if two spins 
are aligned, Bij\ |t) = 0. Likewise for singlet combinations, it holds that Vijd || 
) ~ I it)) = 0- Still, the equivalent of Vij could in principle appear on the nearest 
neighbour bonds. Rotating the spins on one sublattice yields 

VjjVij (afflj - 6|6j)(aja] - bib]). 

The optimal mean field solution — self consistent with lowest energy — turns out to 
satisfy {[aia] — bib])) = for nearest neighbours. The same holds for the equivalent 
of Bij on next-nearest neighbour bonds. Therefore we will neglect these terms. 

The final step is to replace the products of these two operators Vij and Bij by 
their mean field approximations. Set 'jij = \{Bij) for nearest-neighbour bonds and 
Kij = \{T>ij) for the next-nearest neighbour bonds. In the regime of Neel ordering, we 
only consider the translationally invariant solutions, so 7jj — > 7, Kij n. In figure 
|3-1| these mean fields are represented. When the decoupling is made, the detailed 
constraint a^a + b^b = 1 is no longer satisfied as the Hamiltonian will now also map 
out of this subspace. Paradoxically, if the constraint were conserved in detail, the 
mean fields 7 and k would become identical to zero as both alter the number of 
particles on individual sites. The constraint a^a + b'^b = 1 will be relaxed to only hold 
on average, (a'^a + b'^b) = 1. Given the translational invariance of the system, this can 



be replaced by a global version and enforced by a Lagrange multiplier, i.e. 
H Hmf + A ^(a|aj + bjbi - 1). 



(3.4.5) 



Introducing the quantities hp and Ap 
hp = 4:J2K cos px cos Py, 

Ap = 2Ji7(cosp3; + cospy), 
the Fourier-transformed Hamiltonian reads 

n = E, + Y,{hp + A)(aJ,ap + bpb^ - Ap{alblp 



apb-p) 



(3.4.6) 
(3.4.7) 

(3.4.8) 



with 



Er = 2N 



2r 



J2 ( - + ) - A 



A solution to this Hamiltonian can be found through the Bogoliubov transformation 
from fan, 6n) to (ar 



ttp cosh 9p + sinh 9 
ftp sinh 9p + cosh 9 



P' 



P' 



tanh 29r 



Ap 

hry + X' 



(3.4.9) 
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If we define tfie excitation energies Up by 



^{hp + Xy-Al (3.4.10) 

tlie Hamiltonian becomes 

H = E, + Y^ LJpialap + (3p(3l). (3.4.11) 
p 

As in all mean field theories the energy has to be stationary with respect to k, 7 and 
A. This yields the equations 

K = J2 ^2 J" ^ ^^^P"" cos py, (3.4.12) 

7 = ^J2T^(^^^P-^ + ^^^Py)^ (3.4.13) 

1 X — ^ ~l~ A . , 

^ = (3.4.14) 

As only cosines appears in these equations either as single terms or in the right 
combinations, the p = (0, 0) and the p = (vr, tt) terms are completely equivalent. For 
future applications it is also useful to define what will turn out to be the 'condensate' 

m. = (3.4.15) 

This is the sum of the contributions of the p = and the p = (vr, vr) terms in the 
summations ( 3.4. 12| ) and (|3.4.14|) . As these are identical, it is twice the 'p = 0'- 



term. If the system size N becomes infinite, 00, equals twice the 'p = 0'- 

term in (|3.4.13|) . It is no coincidence that also appears as the size of the order 



parameter in equation (|3.3.2|) . By inserting a symmetry breaking term 'riJ2j{~^y^^ 



in the Hamiltonian, it can be shown through a short but subtle calculation that the 
equation ( p.3.2|) — without the staggering (— 1)-' — will hold for the ground state of the 



SBMF Hamiltonian with infinitely small symmetry breaking field rj = 0{1/N). The 
staggering (—1)-' has disappeared from the formulae by the transformation (p.4.2|) . 

As long as there is Neel order the condensate will naturally be positive. It will 
obtain a limit value in the case of an infinite system size N; lim^r^oo^s = const > 0. 
This finishes the discussion of the Neel ordered ground state. 

Next we consider the collinear order. The route to follow is quite similar to the 
one just finished. We can therefore be brief about it. Introduce the transformation 



9+ 
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(a) (b) 

Figure 3-1. The mean fields for the Neel (a) and the collinear order (b). 

Again we perform a rotation over tt around tlie iS^-axis, however this time not on one 
of the sublattices but on every other row. Define the mean fields R = ^(X'j^i+i) over 
nearest-neighbour bonds in the x-direction, 71 = ^{Bi^i+y) over the nearest-neighbour 
bonds in the y-direction and 72 = ^{Bi^i±s:±y) over the next-nearest-neighbour bonds. 
The quantities hp and Ap are now given by 

hp = 2JiKcospx, 

Ap = 2Ji7i cospy + 4J272 C0SP2: cospy. (3.4.16) 



After a Fourier and a Bogoliubov transformation, similar to ( |3.4.8D and ( |3.4.9| ), we 
obtain 

n = E, + J2LUpialap + PpPl). 
p 

This is equivalent to ( |3.4.11| ) where u and are now given by 



iop = J{hp + xy-A 



Ec = 2N^M^l-R'-^) + M^ + 2^l)-Xy (3.4.17) 
The consistency equations for these parameters are 



K = 


1 Zip + A 
— > — cos 

2ujp 


71 = 


1 ^ Ap 

Ar?2^p^°^^- 


72 = 


1 ^ Ap 

/ 7^ COS COS Py 

N ^ 2u)p ^ 


1 = 


1 Zip + A 
p 2cc;p ■ 



(3.4.18) 
(3.4.19) 
(3.4.20) 
(3.4.21) 
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The 'condensate' is defined in a similar manner as before in ( |3.4.15| ), tliat is 

= %±^. (3.4.22) 

The symmetry of tlie collinear order differs from tliat of tlie Neel order. Tliis conden- 
sate contains tlie contribution of the p = and p = (0, tt) terms in the summations 
( |3.4.18|) and ( p. 4. 211) . These are identical and once again the condensate rhs is twice 
the 'p = O'-term . If the system size ^ cxo it also becomes twice the 'p = O'-term 
of (|XT|) and ( gX20D . 

The main weakness of the SBMF approach lies in the handling of the particle 
constraint a^a + 6^6 = 1. It is only conserved on average by use of a Lagrange multi- 
plier, equation ( p.4.5| ). The ground state will have non-zero weight in configurations 
that have either too many or too few bosons on a specific site. This space is different 
from the original spin space and no direct correspondence to the ground state of the 
frustrated Heisenberg model exists. 

A connection with the original, frustrated Heisenberg model can be made that 
involves a restriction of the state and a transformation of the basis; first the ground 
state of the SBMF Hamiltonian has to be restricted to the subspace where the con- 
dition a'l'a + b^b = 1 holds in detail. Afterwards, pairs of bosons [a^a etc.) have to 
be replaced by spin operators. The restricted wave function now lies in the correct 
space and can serve as an approximation to the true ground state of the frustrated 
Heisenberg model. Wei and Tao ||5^ make this connection for the unfrustrated case, 
J2 = 0. The properties that they extract from this approximate ground state agree 
surprisingly well with the numerical results for the true ground state from various 
Monte Carlo calculations. 

In this chapter we will not follow their route back to the spin problem, but remain 
in the larger, bosonic space. The properties we obtain should thus be appreciated as 
such; they are related to those of the ground state of the frustrated Heisenberg model 
and should be seen as indications of the behaviour of the frustrated Heisenberg model. 

The mean field approximation has provided the energies and wave functions of all 
states of the Hamiltonian and we can proceed to the calculation of the spin stiffness 
P.- 



3.5 Spin stiffness in SBMF 

How can ps be defined properly in the mean field Hamiltonian? This question might 
seem trivial, but a closer investigation reveals that there are two approaches: 

• Start with the original Heisenberg Hamiltonian Ti. Induce a twist with wave- 



length q (just like has been done in section p.3|) and afterwards apply the 
appropriate mean field approximation, ps is related to the ground state energy 
of this mean field Hamiltonians, i.e. 
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• Swapping the first two steps of the previous approach; first apply mean field 
and then twist, i.e. 

n /imf /tMFiqj ,Ps — • 

In the mean field approximation we use here, these two approaches give the same 
result. Still, the correct approach is the first one, as there the energy increase of the 
original Hamiltonian due to the twist is approximated. The second approach replaces 
the original Hamiltonian by a mean field one and starts to investigate the response 
of the mean field Hamiltonian to a twist. 

The calculation of the energy (?^(q)) is a simple repetition of the approach for 
q = 0. This allows a fairly direct derivation of the spin stiffness avoiding second order 
perturbation theory. Still, the aim of this chapter is to obtain the finite-size scaling 

relations for the kinetic term T and the current-current correlation J independently. 
These relations will then be used to extrapolate the data of the DMRG calculations 
to infinitely large system sizes. We will therefore follow the -general- route of second 
order perturbation theory. 

We can use the formulae derived earlier in (p.3.6|) , but simply taking the mean 
field ground state to replace the true ground state in this expression will not suffice; 
the second term contains an inversion which we cannot handle in this form. It is 
necessary to insert the mean field Hamiltonian TYmf- In doing so we extend our 
SBMF approximations to the entire spectrum. 

The SBMF approximation thus has to be performed on the twisted Hamiltonian, 
7i(q), and the first and second derivative are needed to calculate the stiffness ps in 
second order perturbation theory. This in itself is very similar to the description above 
for the untwisted situation and left to appendix ^ The expressions for the current 
and kinetic operators are 



F AF 

t = — JijUijivi — rj)(rj — rj){Dlj + — 

^ F 



2 

^ AF 



From this point on, we set q = g(cos0, sin0). Of the two terms for ps in ( p.3.6|) , T is 
evaluated more easily: 

^ = -^(01^0) 

= T7 E J^J^i (r* - rj ) (r» - ) - ^ E '^v 4 " ) (r* - ) • 

AF F 
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For the two orderings considered in the last section, this expression boils down to 

T°^^i = Jl72-2J2/^^ (3.5.1) 
T'^"" = 2 J272^ + Ji (72 sin^ - ^2 cos^ 0) . (3.5.2) 

These simple equations hold for all system sizes A^. 

The derivation of J requires somewhat more effort. First the matrix element 
(0|j|a) has to be calculated. As \a) is an excitation of T^mf, it has to fulfill the 
relation 

|a) = «I-...-aM-...-/3j|0). 

Of these the only relevant ones are 

i«)=«iAio), 

as can be established by applying the Bogoliubov transformations ( p.4.9| ) to Cij and 
J^ij. The matrix element itself is a combination of (0|jFjj|a) and {Q\Cij\a). By means 
of simple algebra, we obtain 

+ fiij\a) = -25p^,p2r5sinh26'p^sin(pi ■ r^), 

i,j=i+rs 

{0\Clj + Cij\a) = -25p^,p2r5Cosh26'p^ sin(pi ■ r^). 

i,j=i+rs 

The Kronecker-delta forces pi = p2 = p, thus Ea = Eq + Up-^ + = Eq + 2ujp. 
Inserting all these expressions in the equation for J, ( p.3.7| ), will give the explicit 
formula for J. The Hamiltonian in that expression is T^mf- For the two types of 
ordering J is: 

I — 2J2nAp cospx)"^ (3.5.3) 
Ji^Ap - 2J272(/ip + A) cospy 
i7i(^P + A) + 2J272(/ip + A) cosp^j ) . 



jneel 
jcoU 



1 ^ sin , ^ , , 
p p 
^rJ2-^( COS sin 



sm (p sm Py 



In the introduction of the SBMF we already defined the p = term of the summations 
separately. The condensate is defined in expression ( p.4.15| ) for the Neel order and 



( |3.4.22|) for the collinear order. Here we also have to be careful with these terms. For 
the infinitely large lattice these equations can be simplified by replacing summations 
by integrals, i.e.: 



^?---^(2^/'P--- 
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but all divergent terms have to be taken separately, e.g. 

1 

N > 1 = + - 



1 /?-p + A 



(2vr) 



dp- 



where is defined in equation ( p.4.15|) . Partial integration over p yields: 



neel 



J' 



coll 



K] COS 



Ji7(ms - 7) - 2J2K{ms 
2J272(m^ - 72) 

+Ji (71 (m^ - 7i) sin^ - ^(ms 

These expressions for J only hold for infinite system size N, in contrast with the 
expressions for T in ( |3.5.1| ) and ( ^.5.2| ). 
The spin stiffness is given by ps 

neel 



T + J, thus for ^ 00 



Ps 
Ps 



_ ( 2J272 - Jl^^ 

2J272 + -^i7i 



m 



(3.5.4) 
(3.5.5) 



In figure |3-2| the numerical results of these formula are presented. The phase transition 
is first order and no intermediate phase exists. This is in contrast with statements 
in the literature |TB|, that suggest an intermediate phase of dimer or plaquette 



order. Einarsson and Schulz studied the spin stiffness on small clusters and ex- 
trapolated those results to the two-dimensional geometry. In their results the spin 
stiffness vanishes in the region 0.4 ^ J2/J1 ~ 0.6. 

Ivanov and Ivanov |^ have applied a different method to obtain p^. They consider 
the correlation function = {Si ■ Sj). Comparison of this correlation function ^ 
with the non-linear sigma model where ^ ~ exp{27rps/T) (here T stands for the 
temperature), yields 



^neel 



ms{Ji'y - 2J2K) 



pf = 7(2^272 - Ji«:)(2J272 + Jill). 



Our expression for the collinear ordering is different from theirs. This is not very 
surprising as they do not take anisotropy explicitly into account. Their result is the 
geometric mean of the two components in ( |3.5.5| ). 

This section was started with a discussion on the route to be taken; whether first 
to twist or apply mean field. It was stated there that both would lead to the correct 
answer. This can now easily be seen. As (J^ij) = (Cij) = there is indeed no difference. 



3.6 Isotropic Scaling 

Often it is impossible to obtain the values of observables in a system of infinite size 
(A^ — > 00). One approach to overcome this obstacle is to derive the values for various 
system sizes A^ and afterwards extrapolate to the infinite size. This is for instance 
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Figure 3-2. The spin stiffness ps in units of Ji as function of the ratio J2I J\ (sohd Unes). 
For the cohinear ordering the spin stiffness in the direction of the ferromagnetic order 
(lower sohd curve) and in the direction of the antiferromagnetic order (upper sohd curve) 
are drawn. The dotted hue is the result found by Ivanov and Ivanov |^6| for the collinear 
ordering. 



done with data from quantum Monte Carlo calculations and exact diagonalisation 
methods. It is necessary to know the size dependence of the observables to obtain a 
good approximate for their limiting values. 



Recently some controversy has arisen about the size dependence of ps |T7|. It was 



suggested that the lack of proper scaling behaviour on small systems for intermediate 
range of frustration 0.4 < J2I J\ ~ 0.6 is an indication of the absence of long-range 
magnetic order. It is therefore worthwhile to take a closer look at the scaling behaviour 
of the various quantities. This is not a hard task as all formulas in the last section 
have an explicit size dependence. We will show that this lack of scaling behaviour 
does not only appear for intermediate range of frustration, but is a general feature of 
small systems. 

We will focus on the Neel ordering as an example. We want to know the scaling 
behaviour of the condensate nisi ci-nd the two terms J and T that make up Ps {ps = 
T + J). The latter two will turn out to have different scaling behaviour. The following 
discussion is entirely based on the fact that the dispersion relation iv^ is a periodic 
function that is smooth and positive everywhere except at and (tt, tt) where we see 
a linear behaviour; e.g. c^k = c|k| for |k| ^ 1 , the spin wave velocity. 



First the condensate is considered. It is defined by ms{N) = (/iq + X)/{Nujq) 



6Spin stiffness and finite-size scafing of the frustrated Heisenberg Model 



or 



N 



(3.6.1) 



In the hmit N —>■ oo the parameters (k^v, Jn, Xn) wiU obtain their limiting value (k, 
7, A) and the summation can be replaced by an integration. Both changes will give 
rise to corrections. We carefully investigate the size of these corrections below. 

It is known that \imN^oofns{N) = iris > 0. This can also be expressed as limTv^oo 
{ho + X)/{Nuo) = nis or inverting this relation into one for uj-p at |p| -C 1 



Ar2 



The two constants in this formula are given by 



1 
2 



{AJ2K,N + Xn) 



X 



N 



(4J2Kjv) 



(4Ji7A^)' 



' AJ2KN + A 



N 



ms{N) 



(3.6.2) 



(3.6.3) 



The suggestive notation (? anticipates that this is the spin wave velocity since an 
antiferromagnet has a linear dispersion relation ujp = c|p| for low energy. For finite 
system size N the smallest q- vector in the summation ( |3.6.1| ) has length |p| = ^ 
(L^ = A^). This means that the q-independent term, Kq/N"^, is small compared to 
the q-dependent term and gives rise to corrections of at least the order (9(1/A^^) in 
the summation. We replace kn, In and Aat by their limiting values k, 7 and A and 
thereby neglect these corrections of order (9(1/A^^). As the term p = is excluded 
the summation still is finite. 

The other effect, replacement of the summation by the integration, can be treated 
quantitatively due to a lemma by Neuberger and Ziman |3^ . They consider a function 
/(p) that is periodic on the Brillouin zone. If this function / (p) satisfies ^ 1 as 

IpI - 



that is periodic on the Brillouin zone. If this function /(p) satisfies 
* and /(p) is non-zero and smooth in the rest of the Brillouin zone, then 

^ l_ f J_ _ ajLyW) 



We respect their notation by taking the argument L'^/W'^. The numerical value a(l) 
0.6208 was computed in []32| . 

In our situation the limit |p| ^ is given by 



1 2cup 

lim — - 

p^o IpI hp + a 



V2x 



'A - 4J2K 



A + 4J2/t 

The quoted lemma therefore leads to 

0.6208 ' 



ms{N) - rris 



N 



IX + 4J2K 

X - 4:J2K' 
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Figure 3-3. The numerical scaling behaviour of rris (soUd hues) for sizes N = 
10^, 20^, . . . , 90^ (bottom-up) compared with the theoretical curve (dotted line) At J2/J1 = 
0.62 the ground state becomes instable and the discussion on the scaling behaviour is no 
longer applicable. 



This result is in excellent agreement with the numerical values we get when the 
equations (|3.4.12| - |3.4.14]) are solved and the obtained k^, 77V and A^v are inserted in 
( |3.6.1| ). This is depicted in figure p-3| . For values of J2 > 0.62 we know from figure |3^ 



that the phase becomes instable and it is no surprise that the data does not longer 
fit the scaling relation. 

The derivation of the scaling behaviour of J proceeds in the same manner as 
above. The starting point is (|3.5.3| ). The result is 



_ 0.6208 V - (4 J2/«)2 0.6208 c 
^/N AV2 VN 4 



Again we find excellent agreement with the numerical results in figure |3-4| upto J2 = 
0.62, where the Neel phase becomes instable. 

To obtain the scaling behaviour of T a more involved reasoning is required. As is 
seen in ( |3.5.1| ) it depends both on kn and •jn- They are part of the set {hn, In, ^n) of 



solutions to ( |3.4.12[ - |3.4.14|) . As these are mutually dependent they have to be solved 
simultaneously. In the formulation of ( |3.4.12| - |3.4.14| ) the divergent 'p = O'-terms 



are still included. We will rearrange these equations to remove the poles. Give the 
equations ( 3.4.12 - 3.4.14] ) the numbers /, // and /// respectively. For the derivation 



of the scaling behaviour we will use the combinations I — III, II — III and AJ2Kj^iI — 
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Figure 3-4. The numerical scaling behaviour of J (solid lines) for sizes N 
10^,20^,40^, 100^ (bottom- up) compared with the theoretical curve (dotted line). 



4Ji7^// + XnIII, or (using ( pX6D , ( pX7| ) and ( pXTOD ) 

UN = -L + T7 2^ ^ (cos cospy - 1) , 

i\ p ZiOp 

,.1^1 

In 
Xn 



1 hp + Xn 
2ujp 

1 + ^ II ^ (^(cosp^ + cospy) - {hp + Xn] 



4Ji7^ - 4J2K^ 



1 



2N 



(3.6.4) 
(3.6.5) 
(3.6.6) 



It is not possible to extract the scaling behaviour from the equations completely in 
analytic form: the exponent of can be found, but the prefactor cannot be derived 
easily. As an example we consider the equation ( |3.6.6| ). 

First we want to replace the set {kn, ■Jn, Xn) in the summation over Up by (k, 7, A). 
From ( p.6.2| ) we know that Up has the general form 



Ar2 



+ ^7V(P)- 



where gN{p) ~c^|ppif-^ < |p| ^ 1. On all lattice points except the origin gN{p) ^ 
Ko/N^. We have uo = 0{1/N) or l/2NuJo = 0{1/N^). For p 7^ we can expand cUp 
in Ko/N'; 



UJr 



9n[P) 



1 + 



Kn 



2N^9Nip) 



+ ... 



9n{p) + 0{ 
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The corrections in the total summation thus become 



p p 



This mean that replacing the set {k^, ■jn, Xn) by (k, 7, A) leads to errors of the order 
0(l/Ar3/2). 

As in previous cases the summation has to be replaced by an integral. Neuberger 
and Ziman provide a lemma well suited for this type of summation. It reads: 



consider the same /(p) as before, 



The numerical value = —0.7186 was computed in Application of this lemma 
again leads to corrections of the order 0{1/N^^'^). In a similar manner the other 
equations for k and 7 can be treated. 

We now have two sources of finite-size corrections. On one hand the dependence 
of K, 7 and A on the systems size. We know the exponent of this correction, 0{N^^'^), 
but not the prefactor. On the other hand we have the correction arising from the 
summations. There we not only know the exponent, 0{N^^'^), but also the prefactor 
in equation (|3.6.7|) . Since both size dependencies are 0{N^^'^) and can therefore not be 
separated, it is not easy to obtain the overall prefactor. However if the size dependence 
of the parameters (KAr,7Ar, Xn) is neglected the following prefactor is found for T: 



^ 0.7186 1 X + AJ2K , , , , 



As can been seen in figure ^-5| this gives a reasonable description for not too large 
J2/J1. The numerical calculations show that at J2/J1 = 0.62 the scaling behaviour 
is of order 0{1/N). In figure |3-5| one therefore observes crossing-over behaviour from 
0(l/A^^/2)-scaling to C(l/iV)-scaling around this point. 

The size dependence of both T and J are now known. The scaling behaviour of 
the spin stiffness ps (= T + J) is dominated by the order C(l/\/iV) -behaviour of J. 



3.7 Scaling for a highly anisotropic geometry 

The DMRG method that we will employ in the next chapter, allows us to study sys- 
tems of fairly large length L while the width W has to remain small. This fact can 
be seen as an invitation to apply a two-step scaling procedure; first the length depen- 
dence is removed from the observables by taking the limit L ^ 00, and afterwards 
one takes the limit W ^ 00. For both steps knowledge of the scaling behaviour is 
essential. The first step is fairly trivial and will be discussed in due time. For the sec- 
ond one, scaling the width W away, some guidance from the SBMF results is useful. 
Both the kinetic term T and the energy density Eq/N can be treated in the same 



66pin stiffness and finite-size scafing of the frustrated Heisenberg Model 




Figure 3-5. The numerical scaling behaviour of T for square systems of size = 
10^, 20^, . . . , 90^. The A^^/^-scale is derived from analytical arguments in the text. The 
dotted line gives the expected behaviour if the size dependence of k, 7 and A inside the 
summations (|3.6.4l3.6.6|) is neglected. 



fashion. The energy density can be simply obtained by evaluating the Hamiltonian 
Ti in the SBMF ground state; the expectation values of the operators can be written 
in the mean fields n and 7. The energy density is given by 

§ = -2J,(2f -1) + J,(2.^-^). 

As we have not spent any time so far on the energy density, we will start with the 
scaling behaviour observed in the literature. Neuberger and Ziman [R^] derive 



N N-.00 N A^3/2^v / 1 

The function {3 is the same as defined in (|3.6.7] ). We thus seek the behaviour of (3 
as L — >■ 00. Neuberger and Ziman's approach can be extended to this situation. 
Still, the algebraic manipulations are quite tedious [^] and using 

00 

C(3) = E - = 1-2021, 

„=i ^ 

they lead to 

(5iL'^lW'^) 1.2021 

lim = 

(W)3/2 27rPy3- 
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As expected all length dependence drops out and only a finite width dependence 
remains. When we neglect the dependence of the parameters k, 7 and A on the 
system size as before in ( |3.6.8| ), the expressions for the energy density Eq/N and the 
kinetic term T are 

El _ ,„„ El = (3.7.1) 



T.-T ^ ^^J^^[A-8J..]. (3.7.2) 



In figure |3-6| we have numerically checked the second of these predictions. Just like 
in the isotropic case, the power is correct and the prefactor is reasonable. 

Unfortunately a similar limit L 00 can not be taken in the scaling expression 
for the current-current correlation J. The reason can be found in the underlying 
assumptions; we supposed that the dispersion relation is linear and the correction 



due to the fact the ur^ has a gap A, as presented in (|3.6.2| ) has been neglected. For a 



highly anisotropic system, L — 00, this assumptions is not allowed. We can illustrate 
this by an example: the summation over reciprocal lattice in the long direction of the 
system can be replaced by an integral and the outcome diverges; 

LWitJM Wi-2nJ \^iTkl 



whereas the infinite size expression stays bounded; 
dk-— = finite. 



(27r)2 J |k 

The analytical expression we therefore should derive, involves the gap A. Above we 
mentioned that it is not possible to obtain the prefactor of the scaling behaviour of 
7, K and A. The gap A is directly expressed in these three quantities and its explicit 
size dependence therefore becomes just as elusive. However, we can resort to numerical 



means. In figure |3-7| we establish that the current-current correlation shows a 1/W 
scaling behaviour. 

To conclude, the width W has taken over the position of the square root of the 
system size , \fN , in all scaling relations when the highly anisotropic limit of L — > 00 
is taken. (V^ = \/LW) This is naturally not surprising as the smallest dimension of 
the system is W which should set the length scale. 

In the next chapter a brief discussion will be given of the behaviour of the energy 
gap A. This is the energy difference between the ground state singlet and the first 
excitation triplet. For an infinite large system with spin waves this gap obviously will 
be zero whereas for other types of order it might remain finite. For reasons to be 
outlined later, we will consider the member of this triplet with 5^ = 1. In the SBMF 
approximation the gap A is given by 
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Figure 3-7. The scaling of a system of 'almost' infinite length L = 4000 and finite width 
= 6, 8, 10, 12, 14, 20 (bottom-up). 
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The factor 2 can easily be understood if we rewrite the operator in the excitation 
operators and /5k; 




The ground state contains no excitations so S^\0) = 0. The = 1-space can be 
reached by creating two excitations; iS^aottolO) = 1 -aoaolO). This is not surprising in 
the view of the fact that the number of bosons should be conserved on a site. Applying 
only one creation operator creates an excitation that does not satisfy this condition. 
When two creation operators are applied to the ground state, their combination 
contains terms that do satisfy this condition, for example, it contains terms of the 
form Ojbj. 

From ( |3.6.2| ) and ( p.6.3| ) we know that the gap A = 2ujo = 2a/K^/A^ has a very 
subtle size dependence as k, X and 7 have an 0(l/iV^/^) dependence but Kq/N'^ has 
an Oll/N"^) dependence. Adding to this the high anisotropy, we can only establish 
the finite-size corrections numerically. The best fit with a simple function is a size 
dependence A ~ ~ ^0) as can be seen in figure |3-8| . This is in line with the 

idea that the gap A can be considered proportional to the inverse of the correlation 
length. The smallest length scale in the system is W leading to a 1/W behaviour. 

Unfortunately, the scaling of the gap does not follow the route that previous 
scaling relations have followed. The scaling behaviour of Eq/N, T and J in the highly 
anisotropic limit can be summarised as replacing the system size N in the formulae for 
the isotropic geometry by W^. For the gap A that would lead to a dependence 
whereas we actually find a l/iW — Wq) dependence. 



3.8 Conclusions 

In this chapter we have employed the SBMF approximation to get the approximate 
phase diagram of the frustrated Heisenberg model. We found two phases, Neel and 
collinear. The energies of both ground states suggest a first order phase transition at 
J2/ Ji ^ 0.595. There is no evidence for an intermediate phase. SBMF approximations 
in line with either dimer-like order always yield higher ground state energies. As 
mentioned before there are several articles in the literature where an intermediate 
phase is suggested 



The second half of the chapter was spent on finite-size scaling. The scaling be- 
haviour of various quantities (p^, J, . . .) was derived for both square, L = W, and a 
highly anisotropic geometry, L ^ W. These scaling relations will be used in the next 
chapter to extrapolate the numerical data to infinite system size. 
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Figure 3-8. For almost infinite length, L = 1000, and various widths W we observe that 
the inverse of the gap A is linear proportional to the offsetted width 1/ ^ ^ W — Wq. 



4 Density Matrix Renormalisation 
Group approach to the stiffness 



4.1 Expressions for ps 

In chapter ^ we studied the spin stiffness in a Schwinger-boson mean-field approxi- 
mation. In that case the order has a distinct orientation and it is clear how to twist 
it. If we do not resort to a similar mean-field approximation and consider the frus- 
trated Heisenberg model on a finite-size lattice, the orientation of the ground state is 
not as well defined. This model has a continuous symmetry; homogeneous rotations 
in spin space leave the Hamiltonian invariant. We also know that in a finite system 
the ground state is unique. It therefore has to be rotationally invariant. Only in an 
infinitely large system spontaneous symmetry breaking can occur and we can asso- 
ciate a direction with the order of the ground state. The only other way to break the 
symmetry is by enforcing an orientation. This is for instance the case in a mean-field 
approximation. 

The orientation of the order is thus homogeneously distributed over the sphere; 
every orientation in spin space is just as likely. The twist we apply has a plane 
associate with it. Previously we twisted in the x — y plane. The fraction of the ground 
state oriented along the z-axis will not be affected by this twist. Leaving out that 
fraction, we only twist 2/3 of the order parameter. To compensate for that we follow 
Einarsson and Schulz iTGl and introduce 



T = -T J =-J 

sym 2 ' sym ' 



with T and J given in (|3.3.6| ). The spin stiffness is given by 



Ps -^sj/m ~l~ Jsym- (^•-^•-^) 

A further complication lies in the numerical nature of our approach. Unfortu- 
nately the ground state we calculate, will not be entirely rotational symmetric. The 
reason for this touches on the very nature of spontaneous symmetry breaking. As the 
Hamiltonian is rotational invariant, it only takes a small field to orient the ground 
state. Basis states in line with this field will prevail. In standard mean-field approxi- 
mations we use this by introducing external fields to fix the orientation. Here we go 
a step further and directly meddle with the basis. By definition we start our DMRG 
calculation with an asymmetric and poor basis. At each step the ground state wave 
function will be symmetry broken in the same orientation as the basis. The follow- 
ing basis truncation will again be asymmetric. Even if we were to start off with a 
symmetric basis, numerical errors would break the symmetry eventually. It is very 
difficult to maintain a global symmetry by means of iterative local basis updates. 
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To repair this partial symmetry breaking, we would have to incorporate the ro- 
tational symmetry exactly in the procedure. The DMRG method allows for certain 
symmetries to be conserved as was extensively discussed in chapter A good exam- 
ple is S^, which we also conserve in the present calculations. This is possible because 
to every basis state we can assign a quantum number s^. The for the entire system 
is then the sum of the of the basis states on the individual parts. It would be nice 
if also the total spin S could be conserved, but we will argue that this is not feasible. 
The ground state lies in the S = space. The conditions for this restriction can 
be easily derived: 

S'\^Po) = ^ E«=.,,,.5"'|^o) = 

^ (^ol (5")' iV'o) = |5"|V^o)|' = 
^ 'S-lV'o) = 0. 

These conditions are rephrased to reflect that we work in the basis that conserves S^, 
S^\^,)=S+\^o)=S-\^o) = (4.1.2) 



A direct consequence of fulfilling ( [4.1.2 ) is a rotationally invariant ground state, as 



the global rotations in spin space over an angle r are given by exp{irS°') 

Our approximation |0o) already satisfies the first condition of ( |4.1.2| ), S^\(t>o) = 0, 



in a standard implementation of the DMRG. Conservation of the second and third 
condition would require for each basis state in an individual part A of the system 
\i)A the image S^\i)A- This would scale up the number of basis states tremendously 
and the calculation would become prohibitively large. We will therefore neglect this 
symmetry and evaluate ground state wave functions that are only approximately 
rotational invariant. 

On the other hand, there is no reason why the symmetry should be completely 
broken; for narrow systems, the DMRG is accurate enough to compensate for this 
symmetry breaking tendency. 

In general the final ground state will be somewhere in between a symmetry broken 
state and a rotational invariant state. The expressions for the kinetic term T and 
the current-current correlation J can be symmetrised to overcome this orientational 
problem; 

tall = JiY^{q- {n- rj)fSi- Sj + J2'Y^{q- {n- fj)fSi- Sj, 

jail = iJiXK^ ■ (^^ X 4' + ^<^2X1(^ ■ (^^ ~ ^ 4' 

and define 

Tall = -^{^o\tall\(po) , Jail = J^iMJall^^^^jalM, 



then the stiffness ps is given by 
1 
2 



a=x,y,z 
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Three current-current correlation J^n have to be calculated for ( [4.1.3|) whereas only 
one Jsym for expression ( [4.1.1|) . In the next section we will see that this is substantial 
more involved and we prefer to use formula ( ^.1.11 ) whenever the accuracy permits us 
to. In practice this means that for narrow systems we use the symmetric form. For 
wider systems the general form is necessary. 



4.2 Calculating wave functions 



The expressions for the stiffness ps in ( [4.1.1|) and ( [4.1.3|) , have to be implemented 



numerically. The first ingredient is the ground state. A standard implementation of 
the DMRG results in a good estimate |0o) of the ground state. The kinetic term T 
can be obtained by a simple measurement on this wave function \(j)o). However the 
current-current correlation J needs a more elaborate approach. 

An inversion has to be performed for J. We will prove that we can invert Eq — H. 
within the subspace spanned by the basis states for the various parts. To derive this 
we first take a step back and inspect the method to calculate the ground state. At 
each iteration of the DMRG the state |0o) in the subspace that has a minimal energy 
Eq = (0o|'^|0o)/(0o|0o), is selected. This energy Eq is always larger than the true 
ground state energy and we thus have a variational principle. At every next iteration 
we can improve upon our estimate by simple minimising Eq further starting with the 
-truncated- outcome of the previous iteration. This variational principle is crucial for 
the method as it enables to distinguish the best approximation to the ground state 
from other configurations in the basis. 

We can design a similar variational principle for the inversion. Define g{x), 

g{x) = ]^{x\n - Eq\x) + (4.2.1) 

where = j|0o) and Eq is the best estimate of the ground state energy known at 
that point in the procedure Q. This function has a global minimum at 

k) = l0i) = ;^r^l0j), (4.2.2) 

as the quadratic term, ^{x\l-i—EQ\x) , is positive definite. In the realm of linear algebra 
|0i) is called the correction vector. This function provides us with a variational prin- 
ciple similar to the one we had before. Moreover the minimum of the function within 
a specific subspace is also given by ( [4.2.2| ) where |0i), |0j) and Ti are now restricted 
to that subspace. The inversion within the subspace is thus the best approximation 
we can make for the global minimum. 



^Unfortunately this is not always the latest calculated energy; it tends to fluctuate. This was 



already discussed at the end of section 2.4 
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At every step the subspace changes and we can get closer to the real inverse |0i). 
If is known with high accuracy, 1/{Eq — 'H)\(j)j) can be obtained with similar 
accuracy A nice by-product is that 



9 {(pi) = ^(0o|j ^^ ^^ j|0o); 



(4.2.3) 



Apart from the prefactor this is essentially the expression for J. 

The basis has to be tuned to present these wave functions {\(po), l^j); |0i)} opti- 
mally. The reason for including the first wave function, |0o); may be evident. The 
other two are necessary, since we need expression (|4.2.3| ) accurately. If the basis does 
not properly represent |0i) or g{(pi) = is incorrect. To adjust the basis 

to these wave functions we have to incorporate them in the density matrix. Let us 
briefly outline the reasoning behind that. Define the truncation error P by 



P 



The tilde denotes the projection of the wave function on the truncated basis. The 
truncation error P has to be minimal. A few linear algebra manipulations leads to 
the density matrix 



Pii' 



P% + Pii' + P\i' ■ 



(4.2.4) 



The density matrix is thus the sum of the individual density matrices. As usual the 
most important states correspond to the eigenvectors of this density matrix p with 
the largest eigenvalues. The different density matrices could have different weights, 
but the effect on the accuracy of the spin stiffness is unknown. We set them therefore 
to be equal. 

The remaining issue of the last section has now also been answered; expression 
( [4. 1.31 ) for ps is more elaborate than ( [4.1.1| ) as three instead of one inversions have 
to be performed. Moreover all these extra wave functions (four in total; two extra 
currents and two extra inverses) have to be included in the density matrix along the 



lines of equation ( 4.2.4 ). Naturally the basis will then be less suited for each individual 
wave function and an overall loss of accuracy will follow. 



4.3 Geometry 

The shape of the systems we study is dictated by limitations of the DMRG. Their 
width is fairly restricted ( maximally 8 sites wide ) and periodicity along the length 
of the system is not feasible. Earlier we explained that the spin stiffness can only be 
measured with the expressions we derived if the axis along which we twist, denoted 
by q, is periodic. We used two different arguments for this, that are both valid in 
their own right; first, the perturbation theory for an open system will not give the 
desired state and energy. Secondly, a similarity with superfluidity exists; a superfluid 
cannot flow freely when there is an impenetrable wall in its way. 



4.4 Scaling 
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Figure 4-1. The square and the tilted square lattice of width W and length L. In the square 
lattice the periodicity is along the vertical lattice-axis. In the tilted lattice the periodicity is 
along the (vertical) diagonal of the lattice. Periodic images of the nearest neighbour bounds 
are depicted by dashed lines 



The system thus has to be periodic in the -narrow- width direction and open in 
the length direction; the shape of a cylinder. 

The model itself puts some extra constrains on the wrapping of the lattice around 
the cylinder. In order to frustrate neither the Neel nor the collinear ordering, the 
periodicity of 2 lattice sites has to be satisfied. 

The two lattices depicted in figure |4-1| fulfill both requirements. The width W of 
the lattice is the number of sites one passes going round the cylinder. The length on 
the other hand is the maximum number of sites one encounters while scanning along 
the long direction. 

When the next-nearest-neighbour coupling becomes dominant, these lattices both 
fall apart in two sublattices. For the square lattice, these sublattices have an effec- 
tive width of W/2, whereas the tilted lattice breaks up in sublattices of width W. 
Knowing that the accuracy of the DMRG rapidly decreases with increasing width of 
the system, we expect results of a strongly decreasing accuracy for the tilted lattices 
with increasing frustration J2. 

The direction q in which the stiffness is measured is different for the square and 
for the tilted square lattice. The square lattice allows a measurement along the axis 
corresponding to the width direction of the cylinder. On the tilted lattice the direction 
to be taken is the diagonal of the lattice. In short: 



Square lattice : q 




Tilted square lattice : q 



1 



V2 



4.4 Scaling 



The DMRG accuracy rapidly decreases with increasing system width W. It becomes 
therefore necessary to apply finite-size scaling theory to obtain quantities of the two- 
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dimensional system. Here we implement a two step scaling where we can use the 
discussion on anisotropic scaling in the previous chapter and the scaling analysis of 
the two-dimensional ITF in sections [I.5.1| and |1.5.2| . 



The first step exploits the strength of the DMRG; for a fixed width W it is nu- 
merically not difficult to vary the length L substantially. By doing so, we can extract 
the dependence of various quantities on this length L and remove it. The remaining 
fraction corresponds to a system of length L = cxo. To obtain this one-dimensional 
scaling behaviour we can not employ the SBMF approximation of the previous chap- 
ter. That was based on periodic boundary conditions whereas the systems considered 
here have open boundary conditions in the long direction. Still, the correct exponent 
can easily be derived. 

While the correlation length ^ is finite, the infiuence of open boundaries on both 
ends only extends over this length ^. The corrections to the bulk behaviour of all 
properties Eq, T, J and Ps is then just a surface term; 

Eo{L,W) Eo{oo,W) _ ^(]_ 



N N \L 

T{L,W) -T{oo,W) = <^{j 

J{L,W) - J{oo,W) = <^{j 

Ps{L,W) - ps{^,W) = Ol^^ 

For the second step, scaling in the width direction, we fall back on the expressions 
derived in the previous chapter. As the length has become infinite, L — oo, it is 
no longer relevant whether the corresponding boundary is open or periodic. We can 
refer to the results for the periodic case derive before; expressions ( |3.7.1| ) and ( |3.7.2| ) 



contain the 0{1/W^) scaling behaviour for the energy density Eq/N and the kinetic 
term T respectively. With the help of ( |3.7.1D we can also extract the spin wave velocity 
c. Moreover figure |3-7| demonstrates the 0{1/W) scaling behaviour that SBMF yields 
for the current-current correlation J. 



4.5 Results 

Reliable finite-size scaling requires a substantial number of data points. We have con- 
sidered various system sizes to make -at least- the first step in the scaling procedure. 



L — i> cx), indisputable. For the square geometry (see figure [4-l| widths W = 4, 6, 8 
were considered and for the tilted geometry widths W = 2,4 were studied. In all 
graphs we set Ji = 1. 

We have used both DMRG variants, the original one proposed by White and 
our implementation |]l5l . For this model we confirm the statements made in chapter H; 
the variant of White is much more fiexible. On the other hand, our variant needs 30% 
fewer states for a similar accuracy in the calculation. Furthermore the ground state is 
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Figure 4-2. The energy density Eq/N as function of the inverse length for width W = A 
on the square lattice. 



more symmetric as the translational symmetry is strictly conserved. A relative small 
extra gain can be made by reusing bases; whenever we start a new calculation that 
differs from the last one in the size of the frustration J2, the bases for the various 
parts for the preceding value of J2 can be used. This reduces the number of sweeps 
needed to about three. 



4.5.1 Scaling to L = 00 

The width W of the system is fixed and for various lengths L the properties Eq/N 
T, J and ps are calculated. Usually we set L = 2W, 3W, iW, 5W. In figures 



and |4-4| we have depicted this for square lattices with W = A and two values for 
J2. Many more lengths L are considered here as it is computationally fairly easy to 
achieve enough accuracy for system sizes up to L = 160. The scaling behaviour of 
0{1/L) is clearly confirmed by the graphs. In the figures is the extrapolated values 
for 1/L = are also depicted. The resulting energy density Eq/N and stiffness ps of 



the infinitely long system are collected for various J2 in figures |4-5| and |4-6| . Figures 
and |4-8| contain the equivalent results for the tilted lattice. 



4.5.2 Scaling to = 00 

In the final, two-dimensional system, the orientation of the lattice, square or tilted, 
does no longer matter. The values of all quantities are equal for both orientations at 
system size LxW = 00x00. However, the prefactors for the finite-size corrections do 
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0.400 



0.390 




Figure 4-3. The kinetic term T as function of the inverse length for width W = 4 on the 
square lattice. 




Figure 4-4. The current-current correlation J as function of the inverse length for width 
W = 4: on the square lattice. 
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Figure 4-5. The extrapolated energy density Eq/N for widths W = 4,6,8 on the square 
lattice. For W = '& ratios J2/J1 > 0.8 were not considered. 
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Figure 4-6. The extrapolated spin stiffness ps for widths = 4, 6, 8 on the square lattice. 
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4.6 Other Indicators 
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not have to be the same. We fit the scaling for both orientations with the same offset, 
but with different gradients. In figures [4-9| , [4-10| and |4-11| this is done for J2 = 0. The 
resulting extrapolations are also plotted in the figures. In all these figure we have 
multiplied the width W of the tilted lattices by a factor of 2 so that both the tilted 
and the square lattices fall within the same range of and 1/W. All data can 

then easily be plotted in the same graph. 

In figures 4-12 and [4-13 we plot the energy density Eq/N and spin stiffness Ps 
for a two-dimensional system. The error bars in these graphs are based on fitting the 
data-points for (cxo, W) to the assumed scaling relations; the errors in the data after 
the first scaling, L ^ 00 are neglected. 

For J2 = there is no sign problem and the literature contains excellent results 



with which we compare in table 3^. The known values for Eq/N, ps,T and J do not 
contradict with our estimates, although the differences are up to 6%. 





Sandvik m 


This work 


Eo/N 


- 0.669437(5) 


- 0.666(1) 


Ps 


0.175(2) 


0.165(10) 


T 


0.3347185(3) 


0.330(2) 


J 


- 0.160(2) 


- 0.165(6) 



Table 4-1. The comparison between our results for J2/J1 = and those by Sandvik |43]. 



Figure [4-12| suggest a first order phase transition as the gradient of energy curve 



seems to change drastically around J2/ Ji ~ 0.6. This is the same behaviour as ob- 



served in the SBMF approximation, figure p^ . 

The error bars of the stiffness ps, figure |3-TB| , increase dramatically while sweeping 
past J2/J1 = |. The reason for this is that the kinetic term T and the current- 
current correlations J for the tilted lattices and the square lattices no longer seem 
to have the same limit in the two-dimensional case. We still enforce this and as a 



consequence the error bars increase dramatically. Einarsson and Schulz suggest 
a region 0.4 < J2/J1 < 0.6 where the spin stiffness vanishes, ps = 0. This is not in 
contradiction with our results although we also can not confirm it. 



4.6 Other Indicators 

In the process of calculating the spin stiffness ps, it is easy to generate the spin-spin 
correlations {Si ■ Sj). Although no finite-size scaling was performed, the correlations 
already given a clear hint what to expect in the intermediate range of frustration. 
For these correlation functions there are no restrictions on the boundary conditions 
as was the case for the spin stiffness ps- To achieve highest possible accuracy we set 
the boundary conditions in both directions to be open. In figure [4-14| the correlations 
are depicted for J2/ Ji = 0.5. To be honest, the Hamiltonian has been modified to be 
more conclusive. Let us explain this. 
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Figure 4-9. The extrapolated energy density Eq/N for square lattices W = 4,6,8 and tilted 
lattices W = 2,4. J2 = 0. The widths of the tilted lattices is multiplied by a factor of 2 to 
get both lines for the square and tilted lattices within the same range. The cross x on the 
axis denotes the extrapolated value. 
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Figure 4-10. The kinetic term T extrapolated. J2 = 0. The cross x on the axis denotes the 
extrapolated value. 
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Figure 4-11. The current-current correlation J extrapolated. J2 = O.The cross x on the 
axis denotes the extrapolated value. 
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Figure 4-12. The energy density Eq/N for a two-dimensional system. This curve is the 
collection of extrapolations done as in [4-9| . The error bars are based on fitting the data to 
the scaling relations. 
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Figure 4-13. The spin stiffness ps of a two-dimensional system. Every point is the sum of 
the kinetic term T and the current-current correlation J obtained by extrapolation as in 



4-lC and 4-11 



In the hterature [g^, ^ there are many suggestions for the intermediate phase. 
Two of these would yields such a correlation picture, namely dimer and plaquette 
phases. They distinguish themselves in a very quantum mechanical manner: the dimer 
phase consists of nearest neighbour singlets nicely stacked next to each other on the 
lattice and all aligned in the same direction. The ground state in the plaquette phase 
basically is a direct product of two vertical singlets plus two horizontal singlets on 
each plaquette. 

Figure |4-14| could correspond to a superposition of two discrete orientations of the 
ground state; one oriented in the vertical direction and one in the horizontal direction. 
Inserting a small perturbation in the Hamiltonian, 

^^ = ^E'5.-5,+,, (4.6.1) 

will lift the degeneracy of the dimer orderings, after which only dimers in the length 
direction will remain visible, whereas the plaquette phase would not suffer severely 
from it. The system depicted in |4-14| has this perturbation ([4.6. 1|) included. It provides 
clear evidence in favour of a plaquette phase. 

As a further indicator, the spin gap A was also briefly studied. This is the energy 
difference between the singlet ground state and the first excitation in the = 1 
space. This excitation is a member of the triplet of lowest excitations. Both for the 
Neel and the collinear ordering the gap A should disappear as there exist spin waves. If 
there is a phase between these two, the gap A might open up. The phases suggested 
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Figure 4-14. For 10x8 system that has open boundary conditions in both directions, the 
correlations {Si -Sj) where i and j are nearest neighbours, are depicted. This is for J2/J1 = 
0.5. A perturbation ( 4.6.1 ) is included to distinguish between plaquette and dimer order. 
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in the literature actually all imply a gap A. The advantage of taking a member 
outside the = space is that it can be considered the ground state in its own 
space and calculated in exactly the same fashion as the ground state itself with the 
quantum number = 1. Unfortunately the energies of both the ground state and 
the excitation grow as whereas the gap remains of order Ji, A ~ Ji. Only for 
W = 4,6 we could obtain enough accuracy to scale away the length (an 0{1/L) 
correction). With only these two points a scaling analysis in the width direction is 
impossible for the simple reason that no optimal fitting can be done (which needs at 
least three points). 



4.7 Discussion 

The scaling analysis of the spin stiffness ps we performed, does not give accurate 
results. It neither supports nor contradicts the existence of an intermediate phase. 
The two successive steps of the scaling analysis have very different degree of success 
and we will discuss them separately; 

The first step, scaling L oo, gives reliable values for the properties of an 
infinitely long cylinder. There are two sources of errors in the properties of this 
cylinder: corrections to scaling and systematic error in the DMRG due to a insufficient 
number of states kept. For the square lattices of width W = 4,6 and the tilted lattice 
of width W = 2, the scaling corrections are the dominant source of errors. The least 
square fits estimate the relative error to be of the order 10"'^. For the square lattice 
of width W = 8 and the tilted lattice for width W = 4, the systematic errors of the 
DMRG should also contribute. The peculiar behaviour of the stiffness in the region 
round J2/J1 = 0.5 for the tilted lattice, figure indicates that this source is there 
even determinant for the overall accuracy. 

The second part, scaling — ^ cxd, is by no means as successful as the first. Figures 
4-9| , [4-10| and |4- 1 1| do not justify the scaling behaviour we assume. Most likely the 



scaling behaviour we derived in the last chapter is not valid for these small system 
widths W. There we also found that for Schwinger bosons this scaling behaviour set 
in at the linear dimension L = 10. 

The literature provides accurate results for the unfrustrated Heisenberg model, 
J2 = 0. Sandvik |^ is the latest in a whole line of authors who have performed 
finite-size scaling on square, periodic systems. Their results fit the scaling behaviour 
nicely and the results are of high quality. While the results are not contradicting ours, 
our inaccurate fit is in sharp contrast with theirs. Naturally we consider an unusual 
geometry (infinitely long cylinders), but the origin of this discrepancy is essentially 
not understood. 

Einarsson, Schulz et al. [l^, ^ performed a similar analysis to ours. Instead of 
scaling the length L 00 first, they considered square lattices with periodicity in 
both directions. From sizes 4x4, 2\/5 x 2v^, 4v^ x 4\/2 and 6x6 they inferred the 
properties of the 2D case. They observe reasonable scaling behaviour in line with the 
unfrustrated case. That their systems are much smaller than ours, makes the contrast 
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with our findings even more striking. 

The stiffness thus does not give a definite answer and we switch our attention to 
the correlation functions. If there exists an intermediate phase, the correlations in the 
ground state clearly hint at a plaquette phase. Zhitomirsky and Ueda ||58| suggested 
before that indeed the plaquette phase is favourable to a dimer one. Still, careful study 
of the dependence of the correlation functions give rise to a few other suspicions: the 
plaquette correlations arise far sooner than the stiffness becomes negligible. Perhaps 
a super solid phase exists? Moreover similar behaviour is observed coming from the 
collinear order, although there dimer (and not plaquette) correlations are appearing. 
This even makes rooms for two intermediate phases. The abrupt change of the energy 
in figure |4-12| suggested that between two of these phases a first order phase transition 
exists. If that is the case, it is most likely that it will occur between a plaquette phase 
and a dimer phase (possibly both with long range magnetic order). 

In the next chapter we will study the correlations further to get insight in possible 
intermediate phases. 
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5 Combination of DMRG and 
Fixed-Node Monte Carlo 



5.1 Introduction 

In chapter Q and ^ it became clear that the DMRG can achieve phenomenal accuracy 
for relatively narrow strips, but the accuracy deteriorates once the strips get wider. 
Section p.6.1| provided a connection between this behaviour and perturbation theory. 
Although this relation only yields an upper bound on the number of states needed 
for a certain accuracy, it is clear that the number of off-diagonal matrix elements in 
the Hamiltonian together with their sizes determine the accuracy. If we still want to 
extract ground state information for wider systems, changes or amendments in the 
method have to be made. 

In this chapter we combine it with the Green Function Monte Carlo (GFMC). 
GFMC is based on the notion of projecting out the ground state. Define the operator 



Q = I — eH with e <^ I |n|, The ground state \ipo) can be found starting with a 



state 10) and letting it relax, 

I^Ao) ~ , n£>l , (5.1.1) 

When e is small enough, not the largest but the smallest eigenvector will be projected 
out by this projection. The Hilbert space is very large 0{2^), making it impossible to 
apply these matrix operations in detail. We want to perform a representative sampling 
of these matrix multiplications. With configurations R, built from individual up and 
down spins a =t, | , \R) = \<7i . . . <Jn), we can construct paths R = Rq, . . . , Rn- 
Insertion of complete basis sets J2r \R){R\ between the individual projectors Q in the 
previous equation, can be read as a summation over paths; 



|z^o) ~ E 1^") 



R l-i=l 



Y[{Ri\g\Ri-i) 



(i?o|0). (5.1.2) 



A path R derives its name from the fact that it only contributes to the ground state 
j-^o) if {Ri\G\Ri-i) 7^ for all successive configurations; only specific paths through 
phase space can be followed. The general assumption of a Monte Carlo simulation 
and the GFMC in particular is that accurate properties can still be obtained when 
only a few of these paths R are semi-randomly selected to represent equation ( |5.1.2| ) . 
As the Hilbert space is very large, 0{2^), even the most extensive GFMC simulations 
can be considered to contain only relative few paths. In practice we always generate 
6000 paths. 

We will apply the GFMC to the frustrated Heisenberg Hamiltonian. Successive 
configurations Ri-i and Ri can differ at most in the orientation of two, nearby spins 
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as the Hamiltonian only contains local spin-pair interactions. The transition strength 
{Ri\Q\Ri-i) is set by the Hamiltonian so that locally a good equilibrium is reached. As 
a consequence the local correlation functions are of good quality. The GFMC uses the 
off-diagonal matrix elements of the Hamiltonian in an essential way to systematically 
probe the Hilbert space. Both these aspect, high quality of local correlation functions 
and the intrinsic use of off-diagonal matrix elements, touch on weaknesses of the 
DMRG; first, as mentioned before the off-diagonal matrix elements severely limit the 
accuracy of the DMRG. Second, the correlation functions are biased through the 
sequence in which the sites of the lattice are incorporated in the basis. When the 
sites of a column are added successively to the basis, the correlations between the 
columns are underestimated. 

The GFMC allows for a systematic bias towards specific paths R without infiuenc- 
ing the expectation values. This is done by means of a guiding wave function 
It becomes a measure of importance for a configuration R and thereby of a path. This 
guiding state |0g) embodies both the greatest strength of the method and its main 
weakness: if a large amount of information on the ground state I'^o) is incorporated 
in the guiding state |0g)) the results will improve drastically. On the other hand, 
without a proper guiding wave function no reasonable results can be obtained. 

Even more emphasis is put on a good guiding wave function when handling models 
with frustration or fermions. These are typical cases exhibiting the 'sign-problem'. A 
good guiding state can maybe not solve this sign- problem, but is can suppress 
it to such an extent that it does not influence the extracted ground state properties. 

With the strengths and limitations of both methods in mind, it seems a logical 
solution to combine them; DMRG can make an initial guess |0o) to the ground state 
I -00)- Although this is a systematic approximation, the local correlation functions 
bear a clear signature of the method. They depend on the mapping from the two- 
dimensional system to a one-dimensional chain that is necessary to apply the DMRG 
(site version). The guess |0o) can improve the GFMC in two ways. Most importantly, 
it can serve as a guiding state, 10^) = \(f)o), to reduce the variance and suppress the 
sign-problem. This guiding state is also used to calculate so-called mixed estimators 
for observables. These mixed estimators also strongly improve with a better guiding 
state. The other aspect, where the DMRG state can help, is in the initial state, 
10) = |0o)) but no matter what the quality of this starting point is, eventually the 
ground state will be reached. Especially the quality of the local correlations will 
increase by this stochastic process. Without a DMRG state l^o); the GFMC would 
require another guiding state. In practice these are relatively simple and consequently 
poor approximations to the ground state, that are involved and complex to construct. 

In this chapter we make the connection between DMRG and GFMC by using the 
DMRG ground state as a guiding state, |0g) = |0o)- First we explain the principles 
of GFMC. Afterwards the sign-problem is discussed and a possible cure is described: 
Fixed-Node Monte Carlo (FNMC) and the extension to stochastic reconfiguration. 
With all that in place we make the connection. In fact the only thing we need to 
extract from the DMRG state |0o) is its value for specific configurations R, (i?|0o)- 
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An algorithm will be introduced to obtain this value for an arbitrary configuration. 
Naturally no table with an entry for each possible configuration can be built as it 
would have a size of (9(2^) just like the number of configurations. An extra section 
is spent on curing a common problem of the GFMC by switching from discrete imag- 
inary time intervals, 1 — eTi, to a continuum, exp(— r7^). This makes the method 
also more elegant. Finally, after all these explanatory and introductory sections, the 
computations are presented and the results are discussed. 



5.2 Green Function Monte Carlo 



GFMC has been widely used for at least two decades now ||2^, [Tl|, In math- 



ematics it finds an equivalent in the Markov chain [25| and in a broader physical 
perspective it strongly reminds of diffusion. 

The method will be explained along the lines of the frustrated Heisenberg model, 
where for the moment we simply ignore the sign-problem. Following sections will be 
dedicated to resolving that complication. The frustrated Heisenberg Hamiltonian is 
a collection of spin-pair interactions 

The last term will not alter any of the spins ai, . . . , in a state |-R) when applied 
to it, the first two terms will allow the exchange of an up- and a down-spin. This 
limits the number of states \R') that are connected to \R) strongly; either they are 
identical or in the case that ai ^ aj, they have spins i and j exchanged, = aj and 
a'j = cTj. Applying the Hamiltonian to a configuration reminds of diffusion as it allows 
the up-spins to hop from one site to the other. 

As mentioned in the introduction we want to project out the ground state lipo) 
starting from a — not yet identified — {(j)) by successive applications of ^ = 1 — eH, 
equation (|5.1.1| ). The wave function (Rlipo) cannot be obtained completely because 
of the size of the Hilbert space. For most physical systems it is even arguable whether 
that is desirable. The physical properties are most important and GFMC focuses on 
the determination of these. 

There are two categories of observables X, conserved ones ([A", 7^] = 0) and non- 
conserved ones ([A", ?-^] ^ 0). The conserved observables A" including the Hamiltonian 
itself, can be measured in a fairly simple manner. The guiding state |0g) will be used 
to construct a mixed estimate with exactly the same expectation value as the required 
measurement. 



It is essential that the observable X commutes with the Hamiltonian Ti as after the 
commutations we use \ipo) ^ ^"^^|0g) ~ for 3> 1. 
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The previous relation, ( |5.2.1| ), does not hold if the observable X is not conserved, 
[A', 7i] 7^ 0. We will show how the mixed estimate differs from the required one, and by 
a simple extension reduce this difference. The guiding state can always be considered 
to contain a component along the ground state and a component orthogonal to it, 

|0g) = |V'o)+^|^i) , (V^ol^o) = (^lIV^l) = 1 , (^i|^o) = 0. 

A good guiding state should have a small perpendicular component, 5^1. When 
the commutations in ( |5.2.1| ) are not allowed, the mixed estimate reads 

(^)mixcd = (^ol'^IV'o) + 5(^l|;f l^o)- 

Relating this to the expectation value of the guiding state, (^gI'^I^g), can reduce 
the corrections to order 

= (^olA^l^o) + 5^ ((^olA'l^o) - (^il^lV-i)) + Oif) 

= {^oW,) + 0{5^). (5.2.2) 

To remove this 0{5'^) term completely forward walking schemes are necessary, but 
it is at present unclear whether this can be combined with stochastic reconfigura- 
tion as there the weights are frequently changed. For our purposes only mixed and 
improved mixed estimates, (A')mixed and (A')improved, are required. 

Next is the description of the stochastic nature of the method. In the introduction, 
paths R through phase space were defined, equation ( ^.1.2|) . A selection from all 
possible paths R has to be made stochastically. If the path R is selected with a 
probability P{R) and assigned a weight M{R), the following expectation value has 
to hold: 



{M^G1<I>) = {X{Rn)M{R)) = Y^X{Rn)M{R)P{R), (5.2.3) 

R 

for all X including X = 1. We use here the local expectation value 



X{R) = 



X\R) 



{<Pg\r) ■ 

The mixed estimate can then be obtained by choosing a large number of paths {-R"} 
and calculating 

_EaX{R:)M{R-) 

\ / mixed 



EaMiR-) 

In practice the configurations Ri in a path R are selected successively. The most 
important advantage of this is that a specific configuration Ri-i connects only to 
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relatively few configurations Above it is explained that Ri and Ri-i can differ at 
most in the orientation of two spins for {Ri\Q\Ri_i) ^ 0. 

The starting configuration, Rq, is chosen according to probability distribution 
Pq{Rq). Each configuration Ri afterwards is chosen with probability P{Ri ^ Ri-i) 
giving an overall probability of 

n 

P{R) = l[PiR, ^ Ri.i)PoiRo), 

i=l 

directly in line with the theory of Markov chains. The probabilities Po{R) and P{R ^ 
R') have to be normalised without any negative elements; 

Po{R) = Y,P{R^R') = 1 , Po{R), P{R ^ R!) > 0. 

R R 

In a similar fashion the weight is also successively constructed from a starting weight 
mo(-Ro) and following weight factors m{Ri) combined with 'signs' s{Ri, Ri-i) and 
the weight of the path is finally rescaled with a factor mfin(-Rn); 



M{R) = mRr,{R„ 



Yls{Ri,Ri-i)m{Ri 



u=i 



mo(-Ro 



A first approach would be to let the Green function Q decide; choosing starting 
positions according to their quantum mechanical probability, = 1), 

and expressing no favour for any specific path afterwards: 



^'^^^'-■'^ e!»^0| - '"(^-)=S«I<«I5I^->I- (5.2.4) 



s{Ri, R, 



:\{R\^. 

{Ri\G\Ri~i) 



i-l 



\{R^\g\R,.i)y 

mfir,{Rn) = {4>G\Rn)- 

The equations above also explain why s{Ri, Ri-i) can be named a sign. These com- 
binations satisfy the condition ( p.2.3|) as can easily be verified using 

s{R„R,.i)m{Ri_,)P{Ri ^ = {Ri\g\R,_i). 

In the implementation it is only necessary to store the latest configuration Ri and 
the weight up to that moment. Given the form of the Hamiltonian where up-spins 
make a random walk through the system, the name walker become suitable for this 
latest configuration. The walkers are thus combined with the weights to yield the 
expectation value (|5.2.3|) . 

In practice too many irrelevant paths are selected with these unbiased settings. Far 
better statistics can be achieved using a guiding wave function (_R|0g). Indeed this 
is the same wave function as was used to complete the mixed estimates ( p. 2.1 ). This 
wave function helps us to distinguish important configurations from less important 
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ones and thereby guides the walkers into the relevant parts of the Hilbert space. The 
easiest way to incorporate the guiding state \(t)G) in our calculation is by defining an 
operator X associated to X by 



R,R' 



R) 



As this is a similarity transformation, the projector Q basically remains the same as 
Q. The set of equations ( |5.2.4D can still be used for the stochastic process replacing 
Q by Q. Only the final and initial weights have to be altered, 

mfin(-R) = 1 and mo(-R) = ^rrr^- 

{<p\K) 

If we choose the guiding wave function as starting position, = |</'g), even the initial 
weight can be dropped, mo(-R) = 1. The transition probabilities P{Ri <— Ri-i) are 
now biased towards the most relevant configurations. This only reduces the variance, 
the expectation values are unaltered. 

The algorithm thus far prescribes that after n projections a — mixed — measure- 
ment is made, new walkers are created and the projections restarts. In practice it is 
far more efficient to continue using the same walkers; the existing set is distributed 
according to ((^g|-R)(-R|V'o) while the initial set is distributed according to (i?|0G). 
Relatively few projections have to be performed to do further measurements that are 
both independent of the last ones and representative for the ground state. 

In this process of successive projections the relative weights of the walkers will 
spread exponentially. It becomes unwise to continue the path of certain walkers with 
negligible weight whereas walkers with large weight deserve extra attention. Still one 
does not want to inffuence the expectation values. 

The technique to perform this task is called branching. Just like the stochastic 
process that replaced the projecting, the two requirements here are that the expec- 
tation values are to remain unaltered and the variance is minimised. 

The easiest approach to choose new walkers out of a set of old ones is to 
draw them from a probability distribution 

Pc 

In case the walker a is selected, the weight of the new walker is set to 



\MJ N 



(5.2.5) 



If one does not want to use correction factors p5|], the weights can even be set to 
unity {Ma 1). Despite its elegance the variance of this method is relatively large; 
it can happen that A^ times the same walker is selected. 

It is possible to reduce the variance of this branching process substantially without 
changing the expectation value. Here we introduce a small extension of the method 



5.3 Fixed-Node Monte Carlo 



95 



introduced by Calandra and Sorella ^ to reduce the variance of the branching pro- 
cess. The essential difference with straight selection of N walkers, is that the scope of 
the stochastic process is limited as much as possible. It contains two distinct steps. 
The first step is not stochastic in nature. Rescale the weights 



M„ = N- 



M, 



a 



This weight is truncated to an integer, int(MQ,). For every a there are |int(MQ)| 
new walkers created, each with the configuration of walker a and with the weight 
defined in ( 5.2.5| ). Once this is done for all old walkers, a set of Nq new walkers is 



formed. The integer number Nq is always smaller than N. On simple grounds one 
expects Nq ^ ^N. 

For the remainder of the weight. Ma = Ma — int(MQ,), the second, stochastic step 
is performed. We assign a probability Pa to each old walker. 



P 



\Ma 



a 



Ea'lMa'l 



The probabilities are put next to each other on the interval [0, 1]. In each consecutive 
interval of length 1/{N - Nq) ([0, 1/(A^ - Nq)], [1/{N - Nq), 2/{N - Nq)], . . .) one 
walker is selected by choosing a random number ^ in that interval and establishing 
to which probability interval Pa this number ^ belongs. In this fashion the remaining 
N — Nq walkers are selected giving a total of new walkers. 

This stochastic part is similar to the method of selecting N times one walker 
out of a set of N old walkers which we described before. The essential difference lies 
in the fact that here one walker is selected per interval. The latter method may be 
less elegant than the original proposal, but it reduces the variance of the branching 
drastically. 



5.3 Fixed-Node Monte Carlo 

Unfortunately, the last section does not tell the entire story. The GFMC is severely 
hampered by the so-called 'sign-problem' in models that contain frustration or fermions 
The frustrated Heisenberg model belongs to this class and with it we will exemplify 
the notion of a sign-problem. 

Figure |5-1| depicts the relevant situation. The Hamiltonian contains interactions 
between |(1)), |(2)) and |(3)). We know that the matrix element connecting |(1)) with 
1(2)) is given by: 

mg\il)) = -e{{2)\n\{l)) = -e^. 

In this move the weight M will pick up a minus sign from S(2),(i). The same holds for 
the move |(2)) |(3)). The move from |(3)) back to |(1)) will also induce a minus 
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l(l)> l(2)> l(3)> l(l)> 




Figure 5-1. The possible exchanges that a up-spin can make with neighbouring down-spins. 
1(1)) ^ 1(2)) and |(2)) ^ |(3)) result from nearest neighbour interactions and |(1)) ^ |(3)) 
arises through the next nearest neighbour interaction. This figure is minimalistic in that 
only a fraction of the lattice is drawn and only the exchanging spins are depicted. 



sign as — e {{l)\n\{3)) = -ef . When returned to the original configuration the weight 
of the walker has thus reversed sign. The sign that a walker picks up following this 
loop 1(1)) — i> 1(2)) — >• 1(3)) 1(1)) cannot be removed by basis transformations or 
guiding wave functions. 

The foundation of GFMC is that if more and more paths in phase space are in- 
corporated, the overall weight J^a^^a increases and the average J2a ^aM^/ J2a 
improves. Here this line of reasoning does not hold; an extra path can suppress the 
previous weights. Given the fact that the underlying stochastic process is a Markov 
chain it is easy to show that the average sign, J2a ^a/ \Ma\, will decrease exponen- 
tially in the number n of projections made. Likewise the signal Y^a^aMa/ Y^a^a 
will become very small with respect to the noise. One can only hope that before 
the noises overshadows the measurements, the ground state value has already been 
reached. Under normal conditions, this is rare, but as will be explained in the next 
sections, one can steer the calculation towards such a situation. 

To complete the argument on the sign-problem, two further assessments have to 
be made. First, in the unfrustrated case, J2 = 0, such a loop as described above, 
does not exist. A basis rotation, S = exp{i27[J2x,y{^ + y)Sl^y) removes the signs all 
together from the projector as 

sgs^ = 1 - sns^ = i-eJiY, sSiS^s^ 

1 



From this equation, it is clear that the prefactors of the non-diagonal terms are no 
longer negative. Once e is small enough this also holds for the diagonal terms. It is an 
example of the fairly general approach of a basis transformation to remove the signs 
from the projector Q. GFMC has indeed helped to establish high quality results for 
the unfrustrated Heisenberg model HO . 



Marshall |]29| has proven that after this rotation the exact ground state \ipo) 



S\ijjo) of the system is free of signs, {R\ipo) > 0. Coming up to our second assessment. 



this proof can be extended to the region of small J2 [^H^ for J2 > we just argued 
that a sign-problem existed, thus there can even be a sign-problem when the exact 
ground state is sign-less! 
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One of the first successful attempts to overcome this problem originated in the 
realm of quantum models with continuous degrees of freedom |jTl|. Later it was ex- 



tended by van An and van Leeuwen [|T|, van Bemmel et al. |^ and ten Haaf et al. 



2l[| to lattice models. These methods are called Fixed- Node Monte Carlo (FNMC) as 
both in the continuous version and in the lattice version the Hamiltonian is altered 
by removing the negative projector matrix elements, 

{B!\Q\R) < or {R'\n\R) > for|i?') ^ \R). 

The sign of the out coming wave function {R\Q^\(j)) is fixed and the weights are 
positive definite. In the continuum version this is all there is to it, but in the lattice 
version the projector Q would be altered so strongly that no clear connection to the 
original system remains. An extra potential has to be introduced to compensate for 
the restriction of the hops. In the description of the fixed-node method we follow 



Sorella ||4g] who made a small extension with respect to the proposal by van Bemmel 
et al. §. 

Define a fixed-node Hamiltonian Ti^^^ according to the following rules: if \R') ^ \R), 
then 

{R'\H^'^\R) = {R'\n\R) if {R'\n\R) < 0, 
= --f{R'\n\R) if {R'\n\R) > 0. 

The diagonal element is offset by a sign flip potential, 
{R\V'\R) = E {R'\n\R). 

{R'\n\R)>0 

(Rin'^'lR) = {R\H\R) + {1 + ^){R\V'\R). 

For 7 = — 1 the original Hamiltonian is completely recovered including the sign- 
problem but once 7 > the projector Q^"^ = 1 — eT-L^"^ contains no signs any longer. 
Van Bemmel et al. 0] considered the case 7 = 0. Note that the bar over the sign-flip 
term, V^^, is only cosmetic, as it only appears in the diagonal terms. On the contrary 
we cannot remove the bar over 7i in the definition of V*^^ as here the non-diagonal 
matrix elements are considered. 



It can be proven that this method is variational [^T], i. e. 

((/.IH^'^I^) - (017^10) = (1+ 7) A(0,0g)>O for all 

with A{(j),(f)G) a well-defined, positive function independent of 7. The most impor- 
tant property for this difference A(0, 0g) is that it vanishes at |0) = \4>g)- A direct 
consequence is that if the ideal guiding wave would be used, |0g) = iV'o), the sign-less 
FNMC would yield the ground state properties exactly (even without the improved 
mixed estimator). 

Within the framework of FNMC we might state that we start with the best 
possible approximation {(pc) that can be made prior to the simulation, and let the 
(fixed-node) Hamiltonian improve on that. 
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This approach can be tested on small systems as there we know both the ground 
state wave function lipo) and its energy Eq\iPq) = H\ipo). Experimentally we have 
found for the frustrated Heisenberg model at J2/J1 = 0.5 that a reasonable guiding 
wave with an energy 

(0g|0g) 

will give rise to an outcome of the FNMC simulation with an error AE in the energy 
that is approximately half of the original error, AE ^ ^{Eg- Eq). More often than 
not this will not do. Only when a gap in the energy spectrum exist and the FNMC 
yields an energy below that of the first excitations, it is clear that the state has to 
resemble the true ground state of the system. In the next section an extension to the 
FNMC is introduced to get substantially closer to the ground state. 



5.4 Stochastic Reconfiguration 

Sorella |^ introduced a method that potentially resolves the limitations of the 
FNMC. He named it Green Function Monte Carlo with Stochastic Reconfiguration 
(GFMCSR). The new ingredient is the reconfiguration. The stochastic part refers to 
branching as defined before. It can be interpreted as a sophisticated method to find 
repeatedly a suitable starting point for a straight GFMC with sign-problem. 

It was mentioned before that the sign-problem does not need to be a great obstacle 
if only a good starting wave function could be chosen. The ground state would 
then be reached before the noise component in the weights becomes dominant. In this 
section, three possible extensions are described, starting with a simple combination 
of FNMC and GFMC and finishing with the GFMCSR. 

The simplest solution would be to target the ground state of the fixed-node Hamil- 
tonian 7i^° first through a FNMC and once that has converged, switch to the projector 
Q. Ten Haaf and van Leeuwen [E^ have performed this routine, with 7 = for Q^""^, 



naming it the power method. A large drawback is that after each measurement the 
routine has to be restarted; new starting configurations have to be generated, dis- 
tributed according to the 

It is actually fairly straightforward to avoid the restart. A FNMC can be set up 
with 7 > 0. To each walker two weights are assigned, and Ma- The fixed- node 
weight is updated as prescribed before using the projector Q^"^. The other, normal 
weight Ma is updated as to refiect the normal projector Q: 

m{R)=m'-{R) but ^(^''^) = ^fP^ ^ s'-{R',R) = l. 

The weights Ma correspond to the projection with Q and will suffer from the sign- 
problem. After n projections, when the average sign J2a ^a/ \Ma\ is not too small, 
measurements are made and a new, sign-less starting point is taken by assigning 
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Still is remains unclear whether the weights Ma have converged enough at the time 
of measurement and more sophistication is necessary. 

There is information in the walkers that both schemes above do not use. In reset- 
ting the weights, equation ( ^.4.1|) , a lot of information on the ground state is lost. All 



kinds of correlation functions just obtained in the measurements are not used to im- 
prove the starting point. GFMCSR, introduced by Sorella [Q, provides a systematic 
method to incorporate this information in the new starting weights. 
If we have a set of observables A"* with expectation values 

the new weights M^''™ should reflect these, 

j:aX\R-)Ma ^ ^ ^ j:aX^iR^)M-^- ^ 

Moreover the average sign, J2a M^ / J2a l^"'^'^!) should have increased substantially. 
The solution to this problem is not unique, but a good handle can be found in the 
flxed-node weights M^. Start with the expression 

with the average 

Ml- ' 

completely in line with previous deflnitions. The prefactors j3i are tuned to satisfy 
equation (|5A2D ||. 

This will yield a starting point with exactly the same properties as observed in the 
last measurements. In a longer calculation one can even consider adjusting the weight 
M^™ to reflect the expectation values A* averaged over several measurements. 

With two weights per walker, branching has to be somewhat different than before. 
The branching is performed on basis of the normal weights Ma and afterwards the 
flxed-node weight M^ is adjusted, = \Ma\- Usually the branching is performed 
just after the reconflguration. 




5.5 A guiding wave function from the DMRG 

In the previous sections we have seen that a good guiding wave function is of tremen- 
dous importance for all variants of the GFMC (straight GFMC, FNMC and GFM- 
CSR). Historically this has been the bottleneck of the GFMC [|^; before a calcu- 
lations could be performed, a large amount of research time had to be dedicated to 
the design of a guiding wave function that would be both similar to the ground state 
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iV'o) and easy to handle in the GFMC. This latter property means that the inner 
product between a configuration (ai, . . . , cttv) and the guiding wave function I^g); 
((Ji . . . cr7v|(/)G), can rapidly be calculated. 

A natural candidate for this guiding wave function is the wave function {R\4>o) 
resulting from a DMRG calculation. This will also overcome the bottleneck, as the 
DMRG is based on a systematic approximation scheme applicable to many different 
systems. Still, we have seen in the previous chapter, that especially for larger width 
W > 8 the state |0o) is quite distinct from the true ground state. The DMRG state 
l^o) in general systematically underestimates the correlations along the length of the 
system for relatively wide systems. 

It is the distance between the DMRG state |</)o) and the true ground state \ipo) that 
the GFMC has to bridge. The obstacles we face implementing the DMRG state |0o) 
as a guiding wave function = (-R|0o) are of a technical nature. The remainder 

of this section will therefore be conceptual straightforward but full of details. 

We want to know the value of the wave function {R\(j)Q), but both memory usage 
and computational effort are an issue. For a single walker the configurations that are 
of interest are the configuration R of the walker itself and those nearby configurations 
{R'} connected by the Hamiltonian, {R'\T-C\R) ^ 0. As these can only differ in the 
orientation of at most two spins, an efficient algorithm can find these values (i?'|0o) 
relatively fast once the value {R\(f)o) is known. We will first show how to obtain {R\4>o) 
and afterwards indicate how to use the intermediate results of this last calculation to 
obtain (i?'|0o) rapidly. 

A DMRG calculation provides a state |0o) on bases of both the left and the right 
part of the system. It moreover gives the necessary transformations to construct these 
sets. Once we have the DMRG state |0o), the representation can be tailored to suit 
our purposes. In appendix the technical details are described. The most important 
modification is to switch to the density matrix basis \a)i for the left / sites and \a)i 
for the right — / sites. The properties of our representation can be summarised as 
follows: 

• For every partition 1 < / < A^ we can represent exactly the same state |0o) as 

a 

This provides us with A^ tables of each m values for y^A|^. It can be seen as an 
extension of equation ( p.2.5|) . There we have only stated that given a partition 



I such a representation can be made. Here we add that for all partitions the 
identical state can be representated in this form. 

All basis transformations A^^^^, for the left and -B^cra' the right part are 
known, 
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\o^)i = E^wl^)l«')m- 

aa' 

The state \a) contains the spin on site / or / + 1 respectively. This wiU yield 
matrices A and B of size 2m^ each. The estimate of the size is clearly too large 
as we neglect the fact the both A and B is very sparse. 

The wave function {R\(f)o) can now be evaluate at the configuration (cxi . . . (Tn\4>o) 
by induction. A reduction to inner products {ai . . . ai\a)i and (o";+i . . . (JatIq;); is made 
via 

(cTi . . . aN\4>o) = J2 ■ ■ ■ (^i\<^)i{(^i+i ■ ■ ■ o"v|a)«- 

a 

Each of these two inner products can be derived inductively; e.g. 

(cTi . . . ai\a)i = A^^,((Ti . . . cri^i\a')i^i. (5.5.1) 

a' 

To optimise the algorithm for the inner products {R'\(j)o), the intermediate results, 
(cTi . . . ai\a)i and (aj+i . . . <Jn\c()i for all a and /, are stored in tables. We can use these 
to readily calculate the inner product of a nearby state \(j'q ■ ■ ■ ctjv) = i^h^h '^^h'^h) 
|(Ti . . . (Jat) with the state |0o) > h)- This new configuration is almost identical to 
the old one apart from the exchange of the spins on sites I2 and h; 

Wo - ■ ■ ^'n) = Wl ■ ■ -(^h ■ ■ - (^h ■ ■ ■ 0"Af)- 

To calculate the inner product, the system can be split up; 

a 

The second part of this expression, (cr/j+i . . . crAr|a;)/2, can be found in the tables. 
The first part can rapidly be built starting from the known, listed inner product 
(cTi . . . o"/^_i|Q;)ij_i and iteratively extending this inner product to location I2 using 

Further substantial reductions can be made. The most important one is to reuse 
most of the intermediate results for both (-R|0o) and {(/2'|0o)} when a walker moves 
from configuration i? to a neighbour R'. Once a walker has propagated far enough 
for the next measurement or reconfiguration and the next walker will be addressed, 
all tables are removed. This is unavoidable as the memory usage has to be limited. 

A typical system is of size L x L = N with open boundary conditions in both 
directions. The calculation of the inner products costs about 2m?N operations for 
the partial inner products of the configuration R itself and N[-\fNw? + Arn?)/2 for 
all others. Here we have again neglected that A and B are very sparse. Still the 
calculation duration will scale as N^^'^m?. If the tables are reused, an extra reduction 
factor of 4 is achieved. 

There is one strong restriction in the wave function of the DMRG; when con- 
sidering a part of size I the density matrix will select states that lies in specific 
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classes. All other classes, ranging from = +1/2 to = —1/2 will not appear in 
the wave function. For the Monte Carlo simulation to relax properly to the ground 
state, configurations \R), that are not contained in the guiding wave function, have 
to be assigned a fixed and small value /3; 

{R\<Po) = "(i?|0o)" = /?. 



5.6 Continuum imaginary time limit 

The guiding wave function is not perfect and this can lead to unnecessary large 
fiuctuations in weights of the paths. If a walker visits a configurations \R) with a low 
'probability' 1(0^1-^)1 ^ 1; which neighbours a fairly likely configuration \R'), 

{R'\n\R) ^ and \{<I)g\R')\ > \{<I>g\R)1 
the local estimate of the energy E{R) gets an excessively large value. 

Only for a perfect guiding wave function we could make the replacement {(pal'H = 
Eo{4>g\ and this problem would disappear. It will have consequences for the weight 
factor m{R) as the projector Q contains the Hamiltonian 7i; 

{(I>g\R) 

Naturally the walker will almost certainly leave this configuration \R) for \R') the 
next projection, but the harm has then already been done. To compensate for this 
situation, one would like to send e — > 0. Without modifications this limit leads to the 
necessity of infinity many projections n. 

Trivedi and Ceperley |^ developed an elegant route out of this trouble. Remem- 
ber that for e <^ 1 

(1 - = e-="^. (5.6.1) 

A continuous time variant can be formulated where only the imaginary time r = en is 
a relevant parameter. Let us describe it for a sign free Hamiltonian, like the fixed-node 
Hamiltonian 7^^^ 

If we start in a configuration the probability to remain in it for An steps is 
given by 



An 
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Both the numerator and the denominator of the expression can be simphfied using 

(113); 



{R\G\R) = exp{-E{R\n\R)), 

j:{R'm) = l-e^^^^^ = eM-eE{R)). 
R' \(PG\R) 

Thus the probabihty is given by 

P{R ^ i?)^" = e-x-p{-Ane{{R\H\R) - E{R))). 
A random number ^ is chosen to set this time, 

Ar = Ane = — - — — -— . 

E{R) - {R\n\R) 

During this time Ar the weight is multiplied by a factor 
m{R)^^ = (j2{R'\0\R?j = exp{-ArE{R)). 

After this time Ar a jump to another configuration \R') ^ \R) has to be made 
according to the transition probabilities 

p^^, ^ p) (^'1^1^) (^'1^1^) 



j:r"MR'V\r) j:r"MR"\'h\R)' 

In this new configurations the walker remains for another time interval. Once the 
total imaginary time r has passed, a measurement, reconfiguration or branch can be 
made. 

The GFMSR does not differ much from the above prescription. The fixed-node 
weight follows it exactly. The normal weight picks up an extra factor during the stay 
of the walker at a specific configuration. 



s{R,R)m{R) 



An 



exp{-AT{R\n\R)) 



exp{-AT{R\n + (1 + 7)V^^|/2)) 
exp(Ar(-^(i?) + (1 + ^){R\V''\R))). 



exp(-Ar^(i?)) 



In the hops the factor s{R,R) = {R'\n\R) / {R'\n^''\R) = 1 or - I/7 is picked up in 
the weights. 

This approach replaces the discrete imaginary time with a continuum and resolves 
the complication of extreme weight factors m{R). 
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Figure 5-2. The sequence in which the DMRG includes the sites in the left part of the 
system, (a) represents the ordinary order, (straight) (b) represents a sequence that is more 
in line with the appearance of plaquettes and all members of a plaquette are added to the 
basis successively (meandering). 



5.7 Implementation issues 

In the previous chapter, the spin stiffness ps was studied. It is possible to obtain 
the same stiffness in a GFMC simulation by a trivial extension of the method by 
Pollock and Ceperley but it is as yet not clear whether the same approach can 



be combined with the GFMCSR. Instead we will focus on the correlation functions 
for various frustrations J2 ranging from J2 = to J2 = 1.0 in steps of 0.1. 

The geometry of the systems are set to 10 x 10 with open boundary conditions 
in both directions since the correlation functions do not require periodic boundary 
conditions as the spin stiffness did. There are three clear advantages of these open 
boundary conditions: first, if a dimer or plaquette phase were to appear, the location 
of the dimers or plaquette will be locked by the boundary conditions; the four corners 
will always contains such a object and the rest of the system can then easily be 
filled in. The second advantage is that DMRG obtains the highest accuracy in open 
systems. Finally, the set of inner products (-R'|0o) can be calculated much faster as 
no neighbouring states exist with one of the first spins exchanged with one of the 
last spins. For all other neighbours the tables with inner products (ai . . .(Ji\a) can 
extensively be used. 



The DMRG states are built in two distinct sequences as depicted in figure |5^ . 
Both are based on adding one site at the time to the basis with m = 75 basis states. 
The usual approach, figure |5-2| (a), is to add column after column, which we name the 
straight sequence. For a plaquette order the meandering sequence of figure ^-2| (b) is 
preferable. The individual sites of a plaquette are then added sequentially allowing 
strong correlations between them. The energy -Edmrg is systematically lower for the 
meandering sequence than for the straight sequence, see table |5-1| (Ful explanation of 
this table will follow in the next section). Therefore we use the meandering sequence 
to build the guiding wave function. With increasing frustration, dimer correlations 
appear in the straight sequence and plaquette correlations appear in the meandering 
sequence. 
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The mixed estimates incorporated in the reconfiguration are the nearest- and 
next-nearest-neighbour correlation functions, 

(iSj ■ ^j-i- 2^ ) [jiixed ) {^i ' '5j+y)mixed ; {^i ' '^i+z+j/) mixed ; {^i ' ^i+x—y) mixed- 

The guiding states effectively share two symmetries with the system geometry: reflec- 
tions in the lines y = 5^ and x = 5|. These symmetries are included in the mixed 
estimates reducing their number by approximately a factor of four. No further geo- 
metrical symmetries are included. The reflection through the diagonal of the system 
is excluded as the guiding states do not share this symmetry. Moreover a dimerised 
state distinguishes itself from a plaquette state by the lack of this symmetry. 

After each reconflguration a branch is performed. We use 6000 walkers en set 



7 = 0.5. Table [5^ lists the imaginary time intervals r between reconflgurations. The 
times r are set to let the average sign J2a ^a/ Z^a \ Ma\ decrease from 1 to about 0.8. 
At the starting of a calculation the average sign tends to drop to a very small value. 
At the start of the computation the conflgurations are fairly arbitrary in during the 
flrst time intervals r they will change frequently. As a consequence the average sign 
at the end of one of the initial intervals will be almost zero and it will only gradually 
increase to 0.8 over about 50 measurements. During this 'thermalisation period' none 
of the calculated mixed estimators can be used for the flnal expectation values. These 
are thus removed when the flnal averages are calculated. 

At J2 = there exists a transformation that will remove the sign-problem alto- 
gether. Although we do not perform this transformation, still the problem becomes 
almost signless. The length of the interval r is set such that successive measurements 
are independent. 

Correction factors introduced by Hetherington |]25| , ^ are not implemented. A 
typical simulation with a guiding wave function built with m = 75 states, takes about 
300 hours on a Intel pentium 300 MHz machine. 



5.8 Results 



Figures |5-3| and p-4| together with table contain the results. Let us flrst describe 



the table. 

For all values of J2 that we compared, the GFMCSR with a guiding state that 
was built meandering through the system, resulted in a lower flnal ground state 
energy Eq. The guiding state itself also has a lower energy that the one obtained 
from a straight sequence. (For J2 = 0.7,0.8 this statement does not hold, but there 
the values are close.) This is a clear indication that these GFMCSR calculations are 
biased by the guiding state. Future research must determine whether this dependence 
can be removed. 

The dimerisations (x-dim and y-dim) indicate whether the translational symmetry 
is broken in one of the two directions; 

2 L/2-1 w ^ ^ 2 -^/^ ^ ^ 

X-dim = _ 2^ '^{^2x,y ■ S2x+l,y) - -^^'Y'Y{S2x-l,y ■ S2x,y), 
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y-dim = . _ . Yl Yl (^^<^y ' ^^,-2y+i) - T77T Yl Y ('^^^.zy-i ■ '5^,2y)- 

V / x=l y=l x=l y=l 

The expectation values (S ■ S) are approximated by the improved mixed estimator. 

The abrupt change of these dimerisation indicators from J2 = 0.6 to J2 = 0.7 
already suggest a first order phase transition at that point. This is in agreement with 
the results of the SBMF theory and the numerical results in the last chapter, figure 









Straight 






Meander 




J2 
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-E'dmrg 




x-dim. 


y-dim. 


-Edmrg 




x-dim. 


y-dim. 





0.3 


-61.30 


-62.33(8) 


0.002 


0.001 


-61.84 


-62.54(4) 


0.012 


0.002 


0.1 


0.06 


-57.96 








-58.53 


-59.25(2) 


0.017 


0.003 


0.2 


0.04 


-54.75 


-56.08(11) 


0.003 


0.004 


-55.48 


-56.22(4) 


0.022 


0.004 


0.3 


0.02 


-51.75 


-53.17(4) 


0.005 


0.007 


-52.50 


-53.38(3) 


0.034 


0.006 


0.4 


0.02 


-49.00 


-50.51(8) 


0.009 


0.015 


-49.92 


-50.60(5) 


0.035 


0.004 


0.5 


0.014 


-46.68 


-47.76(6) 


0.009 


0.058 


-47.78 


-48.34(4) 


0.063 


0.021 


0.6 


0.015 


-45.41 








-46.03 


-46.40(3) 


0.073 


0.022 


0.7 


0.015 


-45.67 








-45.60 


-46.00(2) 


0.011 


0.020 


0.8 


0.02 


-49.16 








-49.13 


-49.60(9) 


0.009 


0.001 


0.9 


0.02 


-53.61 








-53.70 


-54.52(2) 


0.007 


-0.006 


1.0 


0.02 


-58.46 


-59.71(9) 


0.010 


-0.006 


-58.64 


-59.80(8) 


0.007 


-0.005 



Table 5-1. For each degree of frustration the imaginary time interval r, the energy of the 
guiding state -Edmrg and the properties of the GFMCSR state are listed. In the text the 
quantities x-dim and y-dim are explained. 



The figures and |5-4| give a more qualitative insight in this behaviour. The 
average correlation strength between nearest neighbours is negative for each of these 
systems, (j^_i)(vt/_i) J2{ij) {<Si-<Sj) < 0. Depicted are the individual correlation strengths 
relative to the average. The solid lines indicate correlations that are more negative 
and thus stronger than the average. The dotted lines show which correlations are less 
negative than the average. They can even be positive. In all cases the correlations 
indicated with the solid lines are also the largest in absolute terms. To approximate 
these correlation functions the improved estimator is used. Especially in figure |5^ 
^ the ladder structure of the guiding state still persists in the final result. Forward 
walking schemes can reduce this tendency of the guiding state further. Given that this 
infiuence weakens in figure |5-4| , there is no need for this extension to get a qualitative 
picture of the behaviour. 

In figure |3 an abrupt change from plaquettes to dimers appears between J2 = 0.6 
and J2 = 0.7. The same alteration is also observed in the DMRG guiding state (not 
depicted) . 



5.8 Results 
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Figure 5-3. The relative correlation strengths on 10 x 10 lattice. All other nearest neighbour 
correlations can be obtained by reflection these picture in the two dashed lines. The DMRG 
guiding state follows the meandering sequence of flgure ^-2| (b). More explanation is given 
in the text. Reading from top left to bottom right, the values for J2 are J2 = 0, . . . , 0.5 in 
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Figure 5-4. The continuation of figure 5-3; the relative correlation strengths on 10 x 10 
lattice. J2 = 0.5, . . . , 1.0 in steps of 0.1. 
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5.9 Discussion and Conclusion 



In this chapter we have combined the DMRG method with the GFMCSR to analyse 
the properties of the frustrated Heisenberg model. Part of the chapter is methodical 
in nature and in the remainder the physical aspects are investigated. 

On the methodical side, we find that the combination of the DMRG and GFMC 
techniques is successful. A guiding state is generated in a systematic manner. Pre- 
viously it was necessary to dedicate a substantial amount of time to construct an 
approximation to the ground state that is both easy to handle and contains the 
relevant physics. Here, we present a relatively easy alternative. 

The present implementation is not yet state of the art computationally. Neither 
have we used the strongest machines currently available nor is the software fully 
optimal. A GFMC simulation can easily be distributed over many processors without 
communication between them. This allows for the cheapest way of upscaling by having 
many computers simulate independently. 

The Monte Carlo methods themselves could also do with further improvements. 
The mixed estimator will not suffice for accurate expectation values. A forward walk- 
ing scheme § has to be developed for higher precision, but it is unclear how to 
combine it with the stochastic reconfiguration that alters the weights frequently. We 
mention that FNMC does not yield high enough accuracy and an extension is nec- 
essary; the GFMCSR is such an extension but many complications still remain: it 
is evident that the final result is biased by the guiding state. This is even the case 
for conserved observables like the energy, which should be accurately sampled by the 
mixed estimator. Furthermore it remains unclear which and how many observables 
should be used for the stochastic reconfiguration. More investigations are necessary 
in order to answer all of these questions. 

We select the meandering sequence for the guiding wave function as that se- 
quence yields the lowest energies -Edmrg for the guiding wave function and Eq for the 
GFMCSR calculation. The correlation functions of the GFMCSR are qualitatively 
the same as those of the DMRG, although the numerical values lie closer together. 
DMRG calculations of the same system with upto m = 512 support the same qualita- 
tive behaviour although they do not reach the low energy of the GFMCSR. Therefore 
there is little doubt that the qualitative behaviour we find, is correct for the 10 x 10 
system. 

The simulations shed some light on the physical properties of the frustrated 
Heisenberg model. A clear change of the spin correlations presents itself at J2 ~ 0.6. 
This suggests a first order phase transition at that location, in line with both the 
SBMF picture and the extrapolated DMRG results for the ground state energy Eq. 
With regard to the spin stiffness (see figure [4-13|) , it is unlikely that a dimer phase 
exists for stronger coupling then J2 ~ 0.6. The coUinear order must already set in, 
maybe with dimer-like tendencies but still conserving long-range order. On the weak 
couphng side, J2 <C 1, plaquette correlations seems to set in quite soon. The overall 
picture that appears is that with increasing next-nearest neighbour coupling pla- 
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quette correlations build up gradually while preserving long-range antiferromagnetic 
order. At about J2 ~ 0.4 the Neel order disappears in a second order phase transition 
and only the local, plaquette order remains. In a first order phase transition at about 
J 2 ~ 0.6 the system changes into collinear order with dimer-like correlations which 
will gradually fade out with increasing J2. 

Singh et al. |^ recently suggested that the intermediate phase consists of colum- 
nar dimer order where the correlations between dimers inside a column are stronger 
than between columns. We do not expect that to be correct, although figures and 



5^ show similar behaviour. This behaviour is triggered by ladder-like correlations 



in the DMRG guiding states. We expect that a more accurate estimator for these 
correlations will suppress this feature, as the improved mixed estimator still contains 
errors of order eq. ( ^.2.2|) . 



Appendices 

A SBMF Approximation for U\q)nU{q). 

In this appendix the derivation of the Schwinger boson mean-field Hamiltonian for 
the 'twisted' case, q 7^ is performed. In itself this is mostly a repetition of the 
untwisted case combined with a number of small details. In the process the first and 



second derivative with respect to q are also obtained. These are used in section ^ 
to calculate ps- 

The Schwinger boson notation is completely equivalent with the spin representa- 
tion, therefore ?^(q) -defined in ( |3.3.4| )- can be expressed as 



This time we have to distinguish between ferro F and anti ferromagnetic bonds AF 
explicitly as they are not in general simple nearest and next-nearest-neighbours, where 
the Viji^q) and Bij{q) are given by 



The unitary operator W(q) has been defined in ( p.3.3| ). In section ^]3| a paragraph was 



dedicated to the correct treatment of the orientation of the ordering. These conditions 
correspond to the mean-fields here, which have to be taken translationally invariant. 
Without loss of generality we can take them real; 

K,,{q) = ^{V,,{q)), (A.2) 

7ii(q) = \{J3ij{q))- 
The mean-field Hamiltonian TYmf now becomes 

^MF(q) = E •Jij'^iji'i) (^li(q) + ^ii(q) - 2Kij(q 

F 

- E J^iiM (4(q) + %-(q) - 27.,(q) 

AF 

+x E(4a. + blh - 1) - E -^^^ + ^ '^^^j (^-3) 

i F ^ AF ^ 



This Hamiltonian is applicable both to the Neel and the coUinear ordering. The spin 
stiffness has to be derived from the ground state energy of T^Mplq) This is very similar 
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to the situation we had in section p73| The expression ( p.3.6| ) for ps is recaptured with 



3— ^MF(q) 



t 



dq; 



q=0 



q=0 



(A.4) 



and H replaced by T^mf; The states \a) appearing in ( |3.3.7| ) no longer correspond 
to the excitations of the full Hamiltonian but to the excitations of the mean-field 
Hamiltonian T^mf- In order to get explicit expressions for j and t, we will perform 



some algebra. We define Tij and Cij and use (|A.l 



d 

dq 



q=0 



q=0 



-{vi - rj){aia] - bib]) 



Tj){aibj 



bittj] 



The derivatives of neither the mean fields Kjj(q), 7jj(q) nor the Lagrange multiplier A 
will appear in either J or T. It is easy to understand that the dependence on the last 
one, A, can be neglected. This Lagrange multiplier is tuned to make X)i a\ai+b\bi — l = 
0, so the entire terms drops from the expectation value. The first two, K,ij{q) and 7jj(q) 
do not appear in J, as symmetry considerations yield 



q=0 



l(£A.(q)) 



'''-AM)) 



q=0 
q=0 



2 i'^ij) 



2 (^ij) 



0, 
0. 



2 Ndq 
l/_d_ 

Move over in the second derivative these mean field cancel out, e.g. 
d^ 



dq^ 

d^ 



4(q) + ^«^(q)-4'^.:,(q)) + 



VUq)+Vi^{q) -2/€,,(q) 



The expectation value of the first part is zero by equation ( [A.2|) . 
Inserting these quantities in ( |A.4|) gives 



J 

— * 
— * 

t 



F AF 
^ F 



2 k. 



2li 



AF 



These expressions for the current j and kinetic term t can be used to calculate Ps as 
is done in section 



B Transforming the DMRG state to the density matrix basis 
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B Transforming the DMRG state to the density 
matrix basis 

In this appendix it will be shown how the same DMRG state |0o) can be represent ated 
in a density matrix basis for all possible partitions of the system in a left and a right 
part. The left part contains the first I sites and the right part the remaining — I 
sites. 

We define a basis on the left part of the system containing the first / spins 

and let be a basis on the remaining A^— / spins. The state |0o) can be represented 

by 



for a certain /q. If we consider the final outcome of a DMRG calculation, the left part 
will contain all but one site, Iq = N — 1. It is possible to represent the same wave 
function exactly for all different partitions 1 < / < iV. The technique is instructive 
and we will elaborate on it: Define the transformations A\^j^, and Bj^j, by 

ai' 
aj' 

With the help of these transformations we will construct the split-up for Iq + 1. All 
other split-ups are then trivial iterations of the same procedure. Insert the transfor- 
mation of \j)ig in equation ( [B.l| ); 



I0o) = E<^§N)/ob)/o= E <^§^iVN)'.»l/)^o+i 



^§N)/ob)/o= E "^'""'^ 

El/^i')'o+ili')/o+i- 



where the definition \(3ji) = J2i,j,a is used. Orthonormalisation of 

will yield a new basis {|^)zo+i} ^ind a new transformation A^^J^,^. Moreover the prefac- 
tors of the wave function 'Pijt'^ can also readily be deduced. It has to be stressed that 
both the basis and the transformation A''°^,^ follow from this procedure. If 

we were to have an expression for the transformation A^°^^ already, it is replaced by 
this new one. The split-up is changed from lo to + 1 and indeed this approach can 
trivially be extended to yield for all partitions 1 < / < A^ the bases {\j)i} and 

the prefactors 0'^-. 

The next ingredient of our recipe is to switch to the density matrix basis. The 
end of section (|2.2|) and specifically equation ( p.2.5|) explain that this can be done by 



simple basis rotations on the left and on the right basis yielding a representation 

for the left and for the right part. These basis rotations give us furthermore: 
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• A simple representation of the wave function |0o) = J2a y K3i\^)i\^)i- 

• Basis transformations A'-^^^, and B^^^^,. 

This is the notation we will use when deriving the value of the wavefunction (i?|0o) 
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Samenvatting 



Dit proefschrift is gebouwd op drie zuilen: 

• Numerieke berekeningen, 

• Quantum fase overgangen, 

• Quantum spin systemen. 

Van deze drie is de eerste, numerieke berekeningen, waarschijnlijk het meest een- 
voudig te doorgronden. Vaak spreekt men ook van numerieke simulatie om aan te 
geven dat men de werkelijkheid wil nabootsen. De dagelijkse weervoorspelling is 
waarschijnlijk het meest bekende voorbeeld van een numerieke simulatie. Het laat ook 
goed zien, dat deze berekeningen een kunst op zich zijn. Sinds de vorige eeuw weten 
we aan welke regels luchtdeeltjes voldoen. Deze regels, ook wel Navier-Stokes vergeli- 
jkingen genaamd, zijn elegant en eenvoudig, maar de enorme diversiteit aan weer- 
somstandigheden geeft al aan dat het gedrag van de deeltjes hiermee niet doorzichtig 
is. Pas recentelijk kunncn we redelijkc voorspellingen maken. Natuurlijk hebben de 
steeds snellere computers daaraan bijgedragen, maar computerkracht alleen is niet 
opgewassen tegen de complexiteit van dit probleem. We zuilen de computer flink 
moeten helpen met ons fysisch inzicht. Zo kunnen we sneeuwval in het weer rond de 
evenaar uitsluiten. Indien dit soort vereenvoudigingen niet worden doorgevoerd, zal 
de computer vccl tc vecl tijd kwijt zijn om tot bekende conclusies te komen. 

Het is ook verstandig om het doel van numerieke berekeningen te analyseren. 
Het weerbericht heeft een voorspellend karakter. Er zijn evenwel meer redenen om 
simulaties uit te voeren. In de natuurkunde wordt vaak gepoogd de werkelijkheid te 
bevatten in een klein aantal regels en elementen. Aan de hand hiervan stelt men een 
model op. Nu is er vaak onduidelijkheid over de gclijkcnis van zo'n model met de 
realiteit en een numerieke simulatie van het model kan daar Tiitkomst brengen. Men 
kan zich bijvoorbeeld afvragen of de kleur van het aardoppervlak van invloed is op 
het weer. In de praktijk blijkt dat zo te zijn, aangezien een donker oppervlak veel 
meer zonlicht absorbeert dan een licht oppervlak. De verwachting is dus al dat een 
donker oppervlak overdag sneller opwarmt. Is het noodzakelijk om nog meer kleuren 
te introduceren, of kimnen we rood en groen afdoen als half donker, half licht? Met 
een numerieke simulatie kan de invloed van die kleuren worden bepaald om vervolgens 
de nuancering wel of niet op te nemen in het model. 

Verder staan numerieke simulaties ons toe allerlei metingen te doen die in de 
echte wereld niet mogelijk zijn. We kunnen bijvoorbeeld de tijd even vooruitspoelen 
om naar het weer van morgen tc kijkcn, maar we kunnen ook de tijd eenvoudigweg 
stilzetten. In dit proefschrift koelen we ons model af naar het absolute nulpunt, - 
273,15 ° Celcius. Dit is experimented niet mogelijk. Toch kan dit extreme afkoelen 
enorm helpen om eigenschappen van de natuur te begrijpen. 
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Samenvattend zijn numerieke berekeningen nuttig vanwege twee redenen: Ten 
eerste, laten ze allerlei metingen toe die experimenteel niet mogelijk zijn en ten tweede 
kunnen ze een grote steun zijn bij het vinden van relevante natuurkundige regels. Ze 
hoeven evenwel niet eenvoudig te zijn. Om een goede berekening te doen, dient men 
inzicht te hebben in de fysische verschijnselen die men wil analyseren en verder moet 
men in staat zijn om de vertaalslag naar een computerprogramma te maken. 

Ook de tweede zuil, quantum fase overgangen, kent een allcdaagse analogie. De 
fase overgang van water naar ijs door een temperatuursverandering laat zien dat 
eigenschappen drastisch kunnen veranderen bij maar een kleine temperatuursvariatie. 
Deze overgang is niet de enige fase overgang in de natuur. Er zijn vele andere fase 
overgangen bekend. Opvallend is dat ze vaak niets met een temperatuursverandering 
te maken hebben; er zijn vele voorbeelden van overgangen onder invloed van bijvoor- 
beeld druk, dichtheid of elektrisch veld. Indien ijs wordt samcngcpcrst, vcrandcrt het 
weer in water. De temperatuur is hierbij niet veranderd. Als tweede voorbeeld beki- 
jken we de schermen van draagbare computers en horloges. Deze bevatten polymeren 
die normaliter ongeordend en transparant zijn. Zo gauw een elektrisch veld wordt 
aangezet, richten zij zich en worden ze ondoorzichtig. Dit soort overgangen blijken 
vaak op vergelijkbarc wijze tc kunnen worden beschreven als temperatuurafhankelijke 
overgangen. Hoe verschillend ze op het eerste gezicht ook mogen zijn; een duidelijk 
universeel gedrag wordt waargenomen. De overeenkomsten zijn veel omvattender dan 
alleen de gelijke naam 'fase overgang' doet vermoeden. In dit proefschrift worden 
twee specifieke gevallen behandeld, die ieder als voorbeeld kunnen dienen voor gehele 
klassen van quantum fase overgangen. 

In de eerste twee hoofdstukken bestuderen we een model waar een extern mag- 
neetveld de fase overgang bewerkstelligt. Het is zelf mogelijk te laten zien dat het 
magneetveld hier precies dezelfde funktie heeft als de temperatuur in een gerelateerd 
model! 

De andere fase overgang die in dit proefschrift besproken wordt is in de praktijk 
moeizaam te bewerkstelligen. Men moet hierbij denken aan de variabele samenstelling 
van een materiaal. Hier kan niet eenvoudig een vlammetje onder gehouden worden 
zoals bij ijs om een overgang tc veroorzaken. Men moet vele preparaten, ieder met 
een net even andere samenstelling, de revue laten passeren om de fase overgang te 
kunnen bestuderen. Slechts een preparaat heeft precies de juiste samenstelling die 
hoort bij het punt van de fase overgang. Toch blijft de fase overgang ook duidelijk 
invloed uitocfcnen op het gedrag van allc andere preparaten. 

De dcrdc zuil, quantum spin systemen, geeft aan dat we fysische modellcn bekijken 
die uit spins zijn samengesteld. De individuele spins doen denken aan magneetjes met 
een noord- en een zuidpool. Daarnaast voldoen ze aan nog een aantal andere regels 
die diep in de wcrcld van de quantummechanica thuis horen. Bijvoorbeeld ccn spin 
kan in een toestand verkeren die niet duidelijk georienteerd is. Pas als men gaat meten 
zal de spin een unieke orientatie uitkiezen. Dczc spins zijn in ons geval keurig naast 
elkaar geplaatst op een vierkant rooster dat lijkt op een barbeque rooster; elk hokje 
herbergt een spin. 
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In het proces om fysische systemen zo ccnvoudig mogelijk te beschrijven, wor- 
den vaak — pseudo — spins gei'ntroduceerd. De feitelijke deeltjes kunnen bijvoorbeeld 
elektronen zijn, maar alleen de regels die sterke gelijkenis hebben met de regels voor 
spins zijn relevant. Quantum spin systemen hebben dus vaak een voorbeeldfunktie 
zoals we al eerder tegenkwamen in de paragrafen over fase overgangen. 

Nu de zuilen zijn geplaatst, kunnen we ons richten op het fronton. De eerste twee- 
en de laatste drie hoofdstukken vormen ieder een duidehjk geheel. Het doel van de 
eerste twee hoofdstukken is om een nieuwe numerieke methode te doorgronden. Deze 
methode heet Dichtheids Matrix Renormahsatie Groep (DMRG). In hoofdstuk drie 
en vier wordt deze methode vervolgens toegepast op het gefrustreerde Heisenberg 
model. Dit blijkt maar een matig succes op te leveren en in hoofdstuk vijf wordt de 
DMRG met Green Functie Monte Carlo simulaties (GFMC) gecombineerd om toch 
de eigenschappen van dit gefrustreerde model te kunnen analyseren. De belangrijkste 
resultaten van dit proefschrift zijn: 

1. Een goed inzicht in de DMRG is verkregen en gerapporteerd. 

2. Door de combinatie van DMRG en GFMC kunnen we GFMC makkelijker 
hanteren voor een hele klasse van problemen. 

3. We hebben inzicht gekregen in de fase overgangen van het gefrustreerde Heisen- 
berg model. 

De eerste twee hoofdstukken zijn dus methodisch van opzet. De DMRG was al 
zeer succesvol voor spin ketens en wij willen haar toepassen op spin roosters. Daar- 
voor hebben we een zeer bekend model uitgekozen wat frequent wordt gebruikt als 
voorbeeld van een quantum fase overgang. Er is dan ook recentelijk veel numeriek 
werk aan verricht waaraan wij houvast hebben bij onze studie. Alleen kleine systeem- 
pjes met maximaal 32 rijen van 8 spins elk kunnen numeriek worden behandeld. Om 
toch inzicht te krijgen in grotere systemen wijden we de tweede helft van hoofdstuk 
een eraan, de specifieke cffcctcn van de kleine afmetingen te bepalen. Die kunnen 
vervolgens worden verwijdcrd om dc eigenschappen van een oneindig groot rooster 
te verkrijgen. Dit heet eindige-grootte schaling en wordt al geruime tijd toegepast op 
dit soort problemen. De resultaten die we op deze wijze verkrijgen, komen overeen 
met de liter atuur. 

Hoofdstuk twee bespreckt dc DMRG methode en bckijkt wat de mogelijkhe- 
den en beperkingen zijn. Het zict crnaar uit dat de DMRG alleen goede numerieke 
kwaliteit kan verkrijgen voor smallc roosters (strips). Eindige-grootte schaling blijft 
dus noodzakelijk om de eigenschappen van grotere roosters te kunnen bepalen. In 
hoofdstuk vijf zal evenwel een mogelijke uitweg worden gepresenteerd. 

Hoofdstuk drie, vier en vijf richten zich op het gefrustreerde Heisenberg model. 
Ook dit model heeft een voorbeeld funktie. Tevens is recentelijk gesuggereerd dat 
de chemische verbinding CaV^Og er goed door beschreven zou worden. Indien het 
model geformuleerd wordt zonder de quantummechanische aspecten van spins, is het 
gedrag volkomen duidelijk. Deze extra quantummechanische regels geven aanleiding 



124 



Samenvatting 



tot nieuwe verschijnselen die maar ten dele begrepen zijn. De funktie van de tem- 
peratuur rond het vriespunt van water wordt hier bekleed door een frustratie maat 
J 2- Het systeem is maximaal gefrustreerd rond J 2 = 0.5. Voor grotere waarde van J 2 
neemt de frustratie weer af omdat de spins zich dan echt anders gaan organiseren. 
Bij weinig frustratie (J2 ~ en J2 ~ 1) begrijpen we goed wat er gebeurt. In beide 
gebieden gedraagt het systeem zich eender als de klassieke variant. Daartussen zit een 
onduidehjke gebied (J2 ~ 0.5). Waarschijnhjk is daar zelfs een fase die geen klassieke 
evenknie heeft. Dit geeft aanleiding tot het veronderstellen van twee quantum fase 
overgangen. Tussen de onbekende fase en de beide zwak gefrustreerde fasen. 

Om hier zicht op te verkrijgen gebruiken we in hoofdstuk vier de spin-stijfheid. 
Dit is een maat voor het gemak waarmee twee spins, ieder op een ander uiteinde 
van het systeem, in verschillende richtingen georienteerd kunnen worden. AUe an- 
dere spins van het systeem zuUen zich aanpassen aan deze beide spins. Indien spins 
weinig rekening met elkaar houden is het relatief makkelijk deze twee spins anders 
te orienteren; de rest reageert toch niet. Het tegenovergestelde is het geval indien 
de spins wel degelijk rekening houden met elkaar. Deze indicator zou duidelijk an- 
dere waarden moeten aannemen voor de drie verschillende fasen. Aangezien wederom 
alleen kleine systeem numeriek behandeld kunnen worden, is eindige-grootte schaling 
noodzakelijk om de eigenschappen van veel grotere systeem te bepalen. De resultaten 
zijn redelijk te noemen ofschoon hieruit niet met zekerheid kan worden geconcludeerd 
dat deze tussen-fase echt bestaat. 

We kunnen 00k bekijken hoe sterk de individuele spins met elkaar rekening houden. 
Dit is gedaan in hoofdstuk vijf voor een systeem met 10 rijen van elk 10 spins. In 
figuur [4-14| op pagina ^ is mooi te zien dat de spins zich per viertal organiseren. 
(dikke lijnen geven een sterk verband aan). Dit resultaat is niet alleen een duideli- 
jke ondersteuning van het bestaan van de derde fase maar geeft tevens inzicht in de 
onderliggende ordening. 

Hoofdstuk vijf heeft evenwel nog meer te bieden. In dat hoofdstuk introduceren 
we nog een andere numerieke methode, de Green Funktie Monte Carlo simulatie 
(GFMC). GFMC omvat systematisch gokwerk om de eigenschappen van een systeem 
te bepalen. Dit verklaart tevens het tweede gedeelte van de naam. Globaal komt het 
erop neer dat alle spins een orientatie krijgen opgelegd. De onderlinge relaties worden 
vervolgens vastgelegd en de relevantie van deze situatie wordt bepaald. Nadat vele 
situaties zijn bekeken, kunnen door combinatie van de relaties met de relevantie de 
eigenschappen van het systeem worden bepaald. DMRG kan hierbij enorm helpen op 
twee manieren. 

Ten eerste kan DMRG bij voorbaat al inzicht geven in welke situaties het meest 
relevant zijn. De GFMC bekijkt vervolgens alleen deze. Dit bespaart buitengewoon 
veel werk en computertijd. 

Het tweede aanknopingspunt is zeer quantummechanisch van aard. Het blijkt dat 
in modellen met concurerende wisselwerking (frustratie) of met elektronen moeilijk 
een gemiddelde waarde is te bepalen door het zo geheten teken-probleem. Dit prob- 
leem kan verminderd worden door zo veel mogelijk informatie over het systeem in de 
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berekening te verwerken. DMRG heeft erg veel informatie te bieden die in de GFMC 
methode kan wordt betrokken. 

Deze combinatie van DMRG met GFMC kan succesvol worden genoemd. Het 
laatste woord is er nog niet over gesproken maar er zijn een heel aantal fysische 
problemen bekend waar zij uitkomst zou kunnen bieden. Deze problemen varieren 
van supergeleiding tot quantum hall effect. 



126 Samenvatting 



Curriculum Vitae 



Op 27 november 1970 ben ik geboren in Gorssel, een dorp in dc buurt van Dcvcntcr. In 
1989 haalde ik mijn gymnasium j3 diploma aan het Geert Groote College te Deventer. 
Daarna heb ik van 1989 tot en met 1995 natuurkunde en wiskunde gestudeerd aan 
de Universiteit Leiden. Tijdens mijn studententijd heb ik nog enkele jaren geroeid 
op wedstrijd- en recreatief niveau. In 1993 studeerde ik af in de wiskunde bij dr. 
J. Hulshof op een scriptie over particle diffcrcntiaalvcrgclijkingen. Daarna verrichtte 
ik een jaar lang experimenteel onderzoek aan dc Universiteit van Oxford. In 1995 
studeerde ik af in de theoretische natuurkunde bij mijn latere promotor prof. dr. J. 
M. J. van Leeuwen. Dit maal omvatte de scriptie een statistische benadering voor het 
Hubbard model. 

Na een kort verblijf bij een management advicsburcau in Berlijn begon ik in juni 
1995 aan mijn promotie. Het onderzoek zoals beschreven in dit proefschrift concen- 
treerde zich op numerieke methoden om fysische eigenschappen van quantum rooster 
modellen te bepalen. Hiervoor heb ik ook drie maanden prof. dr. S. R. White aan de 
Universiteit van California in Irvine bezocht. Uiteraard woonde ik diverse cursussen 
en zomerscholen in het buitenland bij. 

Naast mijn onderzoekswerk, heb ik nog het werkcoUege 'Electromagnetisme' begeleid. 
Verder heb ik in de zomer van 1996 drie maanden voor een bank in Londen gewerkt. 




0.320 ' 

0.000 0.005 0.010 0.015 0.020 

1/W' 



5000 



4000 



3000 



2000 



1000 












10 



20 
W 



30 



40 



